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ABSTRACT 


A preprocessor is designed to extract a set of features 
that enhance natural clustering by removing extraneous 
information. The design removes time shift and scale 
dependence by taking advantage of invariant properties of a 
Fourier transform followed by a Mellin transform. The 
preprocessor is realized using an FFT and a Mellin 
transform with a conventional error correction term. The 
error term proves to be indeterminate, but the error’s 
bound is identified as the envelope for Mellin correction 
terms. Properties of the Mellin transform are employed to 
rodify the signal so that the error correcting is no longer 
required. The resulting algorithms are tested with 
variously scaled inputs for which closed form solutions are 
known. With a verified modification in place, the 
preprocessor produces features that are invariant to 
Sea? bi ne ed scaling, while retaining enough information to 
classify canonic shapes. A method of improving performance 


is introduced. 
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I. INTRODUCTION AND BACKGROUND 


Pattern classification is the assignment of a physical 
object or event to one of several prespecified catagories 
and is the result of an incomplete theory of perception. 
Although many transducers are availaole for converting 
light, sound, temperature, reflected radar signals, etc., 
to electrical signals, the ability of machines to perceive 
cr to recognize their environment remains very limited. In 
the structured world of communcations engineering, signals 
are designed to be detectable and differentiable. A much 
more difficult problem presents itself when sensing an 
environment through a transducer and récognizing or even 
classifying the elements of that environment on the sensed 
characteristics of the transducer’s electrical output. 
Pattern recognition can be considered a complex 
communications problem (for example, attempting to teach a 
machine to decode signals encoded by nature). It is 
possible to alter the transducer’s output to facilitate 
object classification, but determining how to alter that 
output is not a simple task. The main elememts of a 
classification system is shown in figure 1 [1]. The 
transducer senses, actively or passively, a set of 
characteristics belonging to the OMMieet These 
Characteristics can never be a complete description of an 
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A Classification System 
Figure 1 
9 





object, but represent, hopefully, enough information to 
classify the object as belonging to one of a number of 
classes. For instance, temperature is a characteristic of 
the object and a class, but this feature is of little value 
unless it differs in some way from objects belonging to 
other classes, while the set must include characteristics 
that are common among that class. The preprocessor (or 
feature extractor) aims to reduce the data by measuring or 
quantifying certain properties that distinguish the sensed 
object as belonging to one class and not to others. This 
can be done by discerning key features or using the imput 
to generate another set of features optimized by some rule. 
The values of each of these features is then passed to the 
classifier, which evaluates these features to assign the 


object to a class. 


With varied success, machine pattern classification has 
been applied to a large range of problems/disciplines. 
Fields where it is particularly common tnclude optical 
imagery, acoustic signal processing, radiology, radio 
astronomy, and electronic warfare, to name a few. Work in 
many of these fields was reviewed during the development of 
this thesis and the results derived and demonstrated here 
may, in turn, be applicable to the field in general. This 
effort has been directed toward designing a preprocessor to 


produce a set of features that are invariant to information 


1¢ 








maownee tombe superfluous to classification, but that retain 
enough information to classify an object. The object has 
been sensed by a transducer and has been represented as an 
empirically derived, univariate time series. Such a series 
would be the form of data available from range only radar 
Mmeura which is specifically what the preprocessor is 
@esigned to handle. Returning to figure 1, the object has 
an infinite set of characteristics (here portrayed as ana 
infinite series of discrete values). The 
transducer/receiver has collected some characteristics of 
the target objective in the presence of noise. This 
information is relayed as a set of discrete samples (hi) 
from a band limited signal. The preprocessor is designed to 
determine and code revelent information (Hj) for the 
classifier. If this task was done well, classification 
btecomes a trivial problem. On the other hand, if the 
classifier becomes ideal (capable of resolving an infinite 
number of characteristics in noise) the preprocessor design 
begins to look like a wire. The distinction between the 
preprocessor and the classifier is arbitrary from an 
analytical point of view. When designing a classification 
oyover functionally, a difference is enforced. The 
classifier has little concern for how the features are 
developed, but seeks to efficiently use those provided to 


guess the class of the target object. The preprocessor is 


ee 





probler dependent, needing to produce an optimal set of 


features, Fj, from the sensed data hi. 


A. FEATURE EXTRACTION 


For the purposes of this paper, two generic approaches 
to feature extraction are defined. The 1g IES Ee a 
classification approach, was described above. The second, a 
descriptive approach, tries to define the object in terms 
of the objects’ structural features. This system might 
recognize a car, for example, by breaking up the visual 
picture into canonic shapes, and comparing this to 
previously specified canonic class models. The perceived 
Structure of the physical object iS maintained and should 
reflect the structure of the object itself. This approach 
could be robust to temporary changes in the object itself. 
In the car exarple, knowing that at one end of the car the 
trunk can be opened or closed allows the device to take 
this factor into account. Another important advantage to 
descriptive techniques is that the class characteristics 
may be entered or specified without collecting actual 
transducer generated data to train the rachine. 
Unfortunately, the problem of designing a machine to 
analyse a visual scene to produce a structural description 
Maseproved to be quite difficult. Object description from a 


univariate time series is even more difficult, and if the 


ire 





radiation sensed by the transducer is not from the visual 
Spectrum, the task rapidly approaches the oD OS = Ie For 
these reasons the approach taken was the classification 
approach (to reduce the signal to a set of orthogonal 
meeuures that do not uniquely reflect the structure of the 
meeect, but do retain sufficient information to classify 


the object). 


This paper excludes a detailed description of the 
transducer specification. The problem of the classifier 
itself is viewed as one of partitioning the feature space 
(Ei) into regions; one region for each category. Ideally, 
this partitioning should be arranged so that none of the 
decisions are ever wrong. When this cannot be realized, at 
least the probability of error should be minimized or the 
average cost of errors minimized. The problem is one within 
euayistical decision theory. <Xnowledge of the object 
classes (the transducer and the classifier) are required to 
design the preprocessor, which is the topic of this thesis. 
The preprocessor designed and built here generates features 
from a range only radar video signal. These features are 
used by a general Bayesian learning classifier. The 
supervised learning general Bayesian classifier approaches 
the problem by taking a series of incoming sets of features 
lateled as to their class. From the data, an a posteriori 


density is computed. Each successive set of training data 
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is used to refine the densities” statistics. When the 
classifier has been trained on N classes, the features are 
modified to separate the class volumes in an optimal way 
and to reduce the number of features to one less than the 
number of possible classes. The feature vectors of class 
members are clustered about a simplex point and likely 
boundaries are set up allowing classification of the object 
as telonging to one of the classes, or of an unkown class. 
An N simplex is a collection of N points in (N-1) space 
where the distance between any two of the points is equal 
to the distance between any other two. Thus a three class 
problem transformed into a three simplex in a two 
dimensional plane produces clustering of the three class’s 
about the vertices of an equilateral triangle. The simplex 
Poerdinates are the reduced feature vectors, generalized 
from the training data. In a controlled, simple problem the 
classifier works well, but when encountering real problems 
one class’s feature space will intertwine another’s, making 
it much more difficult to obtain separation in a meaningful 
way. The goal of the preprocessor is to present key 
features that determine class for subsequent optimization 


by the classifier. 
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B. FOURIER - MELLIN PREPRCCESSING 


Common to all univariate, time series classification 
problems are several variables that interfere with the 
meeeenition process. Assuming discrete data processing is 
used, these are addressed in the following order; 
windowing, framing, scaling, sampling rate, quantization 
noise, and sufficient information. for real processing, the 
imput waveform is not sampled for all time. It is sampled 
for a period of time. This windowing of the data corrupts 
the resulting spectrum in two ways [2]. First it introduces 
a periodicity (the inverse of the window length) and 
resulting aliasing to the otherwise infinite spectrum, and 
further distorts the spectrum by a convolution with the 
spectrum of the window itself. 30th of these effects will 
color all of the data in the same way and so can be 
accounted for by deterministic methods.’ Framing can be 
considered characteristic of poor synchronization, whereby 
the pattern of concern is not position stable with repect 
to the window as shown in figure 2a. It seems that even in 
human optical recognition, the eye tends to center the 
pervern prior to recognition, unless trained otherwise. 
Scaling is that property whereby the object field may vary 
in scale or aspect angle, in one dimension or several 
Gimensions as in figure 2b. Before the analog signal is 


sampled, it must be fed through a sharp low pass filter, 
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Shifting and Scaling Variation 
Figure 2 
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because no higher frequency noise can be present without 
being folded onto the valid data. Quantization noise, due 
to the requirement to round off each sample to some 
discrete level, is treated the sare as round off error in 
numerical processing [3]. In all of these problems 
discussed to this point, the effect of this processing is 
to mask the actual feature vectors, raking Palle 
classification system less sensitive to valid pattern 
characteristics. In all recognition problems, it is assumed 
that there is sufficient information present for a pattern 
to be detected and identified by the system. This means 
that there is sufficient variability between the classes, 
but sufficient similarity between those patterns of a 
single class to classify each pattern in terms of those 


classes. 


A verifiable goal of this preprocessor is to produce a 
set of features that are invariant to shifting and scaling 
changes. An apvroach, figure 5, has been proposed and used 
[4-10]. The preprocessor consists of putting the sampled 
data successively through two transforms, a fast Fourier 
transform and a discrete Mellin transform. Integral 
transforms, the Fourier and Mellin included, develop 


naturally from the solution of simple problems in potential 
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An FM Preprocessor 
Figure 3 
18 





theory [11-13]. The Fourier integral transform of the 


waveform h(t), 
Hf) =) Ale) KCL eat (1) 


where K=exp[-2mft], is a principle analytical tool in such 
Gamerse fields as linear systems, optics, probvability 
theory, quantum physics, and signal analysis [13]. Its 
purpose in the preprocessor is twofold, but relies on a 
Single characteristic. The magnitude of the Fourter 


transform is invariant to shifting, h(t-a). 
7 -42mfa p : 
F[Act-a)] =e" 4 (t) ep rMee oe a Se 


[ cs | = | e** Hes f (2) 


This characteristic removes the effects of framing 
inaccuracies, and also permits the averaging of successive 
looks or pulses of data to improve feature resolution, but 
removes much of the information about a signal’s structure 


as discussed later. A discrete Fourier transform, 


ice) = 2hc-B a) an= 6, 1, 2, .00, Meal 
Hee) —S Hin 4) mao, 1, 2% ey NO! (3) 
M =! 
2am n/N 
" _— 
Him) = > Alm) €@ (4) 
“MN =0 


has an equivalent identity, but it is only exact for shifts 


of integer values. 
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= | Hem 7d 748 ol em) | (5) 


Shifts of other than integer values result in errors that 
depend not only on the shift, but on qualities of the 


sampled waveform itself, h(t). 


The Mellin transform is an integral transform with the 


kernel, K. 


i a ie f ba) t-te. = 


is the Mellin transform with respect to the complex 
parameter s=r-j2im. Several simple subdstitutions relate 
this to more common analytical tools. Exponentially warping 
t=exp[x}] changes the appearance of the integral to what is 


often called a double sided Laplace transform [14], 


Hess =| £cet) eala (3) 


The transform is invariant to t domain scaling when taken 


with respect to the imaginary part alone, 


sit a Kx) 8704 1g (9) 
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Equation (€) is recognized as the Fourier intergral of an 
exponentially distorted waveform h°(x). The modulus of the 
expression on the right is the magnitude of the Fourier 
meanstorm of the exponentially distorted time function. The 
property to be exploited in a Fouriter-Mellin (FM) 
Peeprocessor is that the modulus of the transform ins, is 
invariant to t-scaling. Given a time waveform h(t), its 
Mellin transform H(s) is given as equation (7). Scaling the 


entire t-domain by k, 


mA e/e)] = ( GVA) tale 
Q 


(12) 
Letting r=t/k, and remembering that s is imaginery, 
Z aS 
ke { AorrittAhr = Be H¢s) 
[Ao Hes) | =)! secsp) | (11) 


Much effort and detail is spent implementing equation (7) 
Mertally in Chapters [i and III. The rest of this chapter 
is devoted to providing some required background on the 
€iscrete wersio2s of the Fourier transform and some 
properties upon which the FM preprocessor depends. The 
treatment here is brief, being mainly a review of basic FFT 
concepts and as such may ve skipped without loss of 


content. 


A discrete Fourier transform is not computationally 


efficient and so leads to impractically long processing 


al 








times. The fast Fourier transform (FFT) efficiently 
caemputes the discrete transform and is used in the 
preprocessors built for this thesis. Other properities of 
the FFT should be presented before getting into the 
detailed design of the preprocessor. The first and last are 
concerned with symmetry. If h(m) is real, as in the case of 
the sampled data, then the frequency spectrum of that data 
is even, !Fe(n)}!=!He(-n)!. H(n) has a real part and an 
imaginary part, Re(n) and Im(n), while h({m) has an even 


he(r)=he(-m), and an odd ho(m)=-ho(-m) part. 


=! 
Re Cm) = Ss he Can) ¢os (28 mom ) 
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The odd part of h(n) times the cosine kernel summation, and 
the even part of h(n) times the sine kernel summation are 
both zero. H(n) can now be seen to have an even and odd 
part. Taking the magnitude of an odd function makes it 
even, proving that the Fourier transform of a real series 
is a spectrum whose magnitude is symmetric about f=0. This 
is of course true for both the integral and discrete 
Fourier transforms. The second line of symmetry is an 
effect of discretizing the signal and its spectrum for a 
discrete transform as evinced by its theoretical 
development. Before froceeding though, the convolution 


theorem is required. 


Ge 





The convolution of the two functions h(t) and g(t) is 
defined as the familiar integral, 


gy Ce) = (her) g (4-0) Aap = £4) tS) (4) 
=f" ge Klt-m dn = ci) £2¢t) (13) 


The relationship between convolution and the Fourier 
integral is very important to modern analysis and 
contributes to making the Fourier transform a key analytic 
tool. The convolution theorem states that if h(t) has a 
Fourier transform F(f) and g(t) has the Fourier transform 


G(f), then h(t)*g(t) has the Fourier transform H(f)G(f). 


y Ct) = ZL Ct) x5 (4) 


2 aD ig Hit) GCF) (14) 


Proving this, the Fourier integral is used directly on the 


convolution integral. 
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From the shifting theorem already presented, 
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And finally, 


rod = F[gCWehcal= Gee) te) (17) 
It can be shown similarly that, 


F'C G4) ¢ He] = 9 ACE (18) 


Clearly, convolution in one domain is simple multiplication 
in the other domain. Although not needed for the pending 
development of the discrete Fourier transform, another 
important relationship, known as the correlation theorem, 
can be appropriately dealt with here. The correlation 
integral, 
2 

3) = J Le gree at (18) 
has an operation with which it forms a Fourier pair as did 
the convolution-multiplication operations. This theorem can 
be established as before, 


4 Rt 
eo? 


LCF) “So [ {hong crsedetr | “At 


= seal hon) cos Canf adder t4 if Wy) cin Cantn dele | 
oe oO (20) 
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The term in brackets is the complex conjugate of H(f) and 


is denoted by B¥(f} in the final form of the theorem below. 


F[ [ Lorgervae] = GCF) H*¢#) | 
(21 


We will now continue on to the discrete Fourier 
transform starting with a continuous waveform h(t). The 
waveform is sampled or multiplied by a string of delta 
functions, s(t). 

CO 
h(t) s Ct) = >, Alma CE ape A) ey 
an 2 & 
where capital delta is the sampling interval. The infinite 
sum is not realized and must be windowed, in this example, 
by w(t)=1 for @6#t #(M-1)4 and zero elsewhere. So that 


now, 
m-!l 
Ri) s(e)wle) = D>” RemA) I(e-md) (23) 


an 0 
Recalling the convolution theorem, the multiplications in 
the time domain corresrond to convolutions in the frequency 
domain with the following results. H(f) is convolved with 
the window functions” spectrum and will have the apparent 
effect of introducing ripples because of the window’s 
Significant sidelobes. The rippling may be minimized by 
icesime duwindow functlonwith small sidelobes at the cost 


of other, perhaps more acceptable, spectral degradation. 
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Hee? is convolved with the sampling function and 
S(f)=I(f-n/4 ) has made the spectrum periodic with respect 
to the interval F=1/4a . The spectrum coming from ae real 
waveform is first symmetric about f=@ as shown before. Now 
because of its periodocity the spectrum is symmetric also 
about f=F/2 (the Nyquist or folding freauency). One final 
step remains. The Fourter spectrum is also taken at 
discrete points 1/T apart. The result in the time domain is 
the convolution of the sampled, windowed signal with 


I(t-nT), which is a periodic signal with T as its period. 


The FFT algorithm used in the preprocessors developed in 
this thesis uses a Cooly-Tukey, base two algorithm [15]. 
This is documented where it is used in programs included in 
the Appendices. The algorithm uses N samples where N is two 
to an integer power. In the transformed domain, due to the 
symmetry shown above, there are N coefficients (only N/2 of 
which are unique as shown in figure 4). The original 
waveform must be band limited prior to samvdling to minimize 
aliasing. The resulting frequency spectrum should approach 
zero at the folding frequency where up to hal? of the power 


can be aliasing noise. 


In the following two chapters, several different 
preprocessing algorithms are developved and their 


performance compared with canonic inputs. In Chapter IY the 


ae 





pemeen Of Ship identification with range only radar is 
discussed briefly and one preprocessor is applied to 
several ship profiles. In the fifth and final chapter, 


conclusions are drawn from the work done here and follow-on 


efforts are recommended. 
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II]. DIGITAL MELLIN TRANSFORM 
BY EXPONENTIAL WARPING 


in the first chapter, the rodulus of the Mellin 
transform was shown to be invariant to scaling. A detailed 
exarination On the mechanics involved suggests an 
implementation that is widely used. The Mellin transform of 
a t-domain function h(t) is given in equations (7) and 


(9°). 


Hsp = f Lee) e* de a 
Hs) -S £cae*) e*"olx Kee ) 


Delta has been added, corresponding to the sample interval 
in the t-domain. Again it is noted that (9°) is a Fourier 
transform, where s=-jet@m. Solving the integral for the 


effect of a t-scaling ty the factor k, 


© /n 09 Ss 
("£cem) tol = \ Ree x bn Rk) 1p 


=e shit Hes) (24) 


Clearly, in the Fourier integral, the scaling factor k has 
become a shift for which the modulus of the transform is 
invariant. The exponential warp alone has transformed the 
scale factor into a shift. A prerenguisite is that the 
t-domain signal has no shift. If there were a shift, it 
does not transform to a simple factor or shift in the 


ac 





Mellin domain and so cannot subsecuently be removed by 


taking the magnitude of a Fourier transforr. 


Mop lementatiwen of a discrete Mellin transform is as 
difficult as it is with the Laplace transform [13]. Once 
transformed, characteristics in the new spectrum are 
adifficult to relate to the original signal characteristics. 
One hypothesis relating the two domains is generated by 
comparison to the Fourier transform. The power spectrum 
associated with the Fourier transform can be used to detect 
Pemmodicities in the physical function, since the wave 
numbers at which sharp peaks of the spectrum occur give the 
wavelength of such periodicities. By analogy, the positions 
of the peaks in the spectrum associated with the Mellin 
transform are said to give the magnification or compression 
Which will produce features in the physical domain. 
Further, this stretching and compressing is identified as 
periodic in nature [21]. This seems unlikely because the 
Mellin is invariant to magnification/compression and does 
not behave well (a scale factor that is a function of t and 
k(t) as seen in (24). The Fourier spectrum models the 
original signal by a set of weighted sinusoids of varying 
frequencies, and therefore, naturally display periodicity 


and is invariant to shift. The Mellin is also a weighted 
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sum of sinusoids, but whose magnitudes are inversely 


weighted by ¢t. 
= a arte mt 
Hos) = | 40 a At (25) 





Values for h(t) for @ t 1 are far wore important to the 
sum than those beyond that point. This characteristic 
contributes to the difficulty of realizing discrete Laplace 
and Mellin transforms and iS a major topic covered in this 


chapter. 


The numerical approximation of the Fourier-Mellin (FM) 
transformation by exponential warping is functionally 
diagrammed in figure 5. A FORTRAN program using the 
algorithm described in this chapter is included in Appendix 
A. Referring to figure 5, the input samples are from a 
pulse whose duration is finite and less than that of the 
Sample window. For this case, no spectral distortion is 
experienced by filling zeros behind the sampled data. The 
only effect of the zero filling is to add spectral 
resolution in the frequency domain. Once through the first 
FFT block, and the magnitude taken, Ef is symmetric about 
zero and at the folding frequency. Thus, ror N time samples 
and filled zeros, there will be N/2 unique spectral 
coefficients. This unique spectrum is interpolated, 
resampled as a warpved signal and fed to the final FFT 


block. A correction is added, and the modulus is taken. The 
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resulting FM features are invariant to shifting and scaling 
ier des time domain. The FYP block was covered in sufficient 
detail in the vreceding chapter. Some effort will be spent 
in discussing the warping itself and the need for, and the 


development of, a zero point correction. 


A. ALGORITHM DEVELOPMENT 


This section uses its own notation to address the 
requisite exponential sarpling. The series to be 
transformed is h(f). Its Mellin transform is E(m) where r 


is the Mellin frequency in s=-jef*’m. The transform is, 


Hoan) = J Ace) tet (2€) 


Letting f=Fexp{x] as before, where F is added corresponding 


to the sample interval 
Hom) pi hk CF e*) eX dby (27) 
The need is to evaluate M equally spaced samples in x, at 
ORG ey CM-I) Xx (28 ) 


while the data consists of N equally spaced samples in f, 


Cre ee to aC NGI) F (29) 
Assuming that F is a small enough interval to properly 


Characterize h(f), it is sometimes advisable [2,12] to 
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choose the exponential sampling interval (X) such that the 
largest intersample spacing in h(Fexp[x]) is equal to kF 
where k=1. The other set of conditions used to uniquely 
specify the new samples are that f=F and x=@ will be the 
lowest limit for interpolation, while (N-1)F and (M-1)X are 
equated as the upper limit. The airs t requirement 


constrains the choice of M by 


CM mt) X = 
/ - e aq ~ 1) esd Pe 2x (3) 


Meeting the second condition, the end points in each 


sampled series are equated yielding, 


Substituting (31) into (3@) while applying the exponential 


series approximation gives, 


PS ee y . x 
e = \ fa aa | eae spe 
va we (32) 





Ne=.f 
One more substitution, (32) back into (31) sets up the 
desired result where M is now specified to exponentially 
Sample from f=F to f=(N-1)F with kF being the largest 


interval between samples. 


M = i A ON-1) (33) 
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As N gets larger, M explodes. If N=1E€ then M=41, if N=32 
then M=1@6, if N=64 then M=261, and so on. This strict 
requirement can be corpromised depending upon the 
application. The other extreme [18] is the criterion that 
requires the smallest interval between samples to equal the 
interval between uniformly spaced Nyquist sampling, with 
the intervals increasing exponentially thereafter. The 
specification permits analyzing frequencies approacning the 
largest which can be analyzed with unifcorm spacing. This 
greatly reduces the number of samples and decreases the 
required procesSing time, but is limited in application. 
Two factors mitigate the stringent requirements imposed by 
(33) where k=1. First, using N unique, uniform samples 
yields N/2 unique, uniform samples in the FFT domain. A 
related consideration is that the values of the spectral 
components approach zero at the folding frequency because 
the original signal has been band limited and some over 
sampling is normally recommended. Secondly, the inverse #f 


weighting apparent in equation (25) 


tn j lean be 
Hom) = { - ea" sist (25°) 
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attaches a decreased importance to h(f) near the folding 
frequency fn. These two effects combine to permit a much 
lower sampling rate once the modulus of the FFT has _ been 


taken. In spite of this, the stiffer rule (33) is used for 
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this work to generate the best 'M dorain possible. Whatever 
rule is used, once the exponential sample points have been 
computed (0). an FM preprocessor, they needn’t be 
Mmeeomputed, but can be stored for rapid access during the 
interpolation. Some interpolation must be performed to 
approximate the spectral values of the new sample points. 
Third or even forth order Lagrange polynomials have been 
recommended and used for this purpose with apparent 
success [9,18,19]. The advantage of the Lagranze technique 
Mmemenaat 10S notation 1S particularly compact, consisting of 
Simple summations and repeated ProOocuc.s that lend 
themselves to digital realization. Unfortunately, for the 
data sets used in this thesis, the third order algorithm#r 
was observed to behave poorly, adding a ripple in regions 
of rapidly changing slope. That this might be the case was 
suggested by the advice that Lagrange is very good near the 
central data point when the order of the polynomial is 
known to be the same as the order of the approximation, 
Otherwise it is test left alone [29,21]. In its place, a 
second order polynomial was used to interpolate the warped 
Samples with results that were nearly indistinguishable 


from the actual waveform, as seen later in figure 8. 
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3. ZERO POINT CORRECTION 


Another problem becomes apparent when the Mellin 


transform of h(f) is recelled, 





Hes) a it si flat { Ace) ex (34) 


where s=-j2mrm., The exponentially sampled waveform 
described above is applied to an FFT block. As f approaches 
the folding freauency, h(f) tends to zero. Unfortunately, 
as f approaches zero, the value of h(f) is not zero. In 
fact it is frequently rather high. To make matters even 
more interesting, the left integral in enuation (34) 
clearly shows that the closer to zero f gets, the more 
important h(f) becomes to the integral. Several solutions 


to this problem are considered below. 


One practical, simple apprcach is to set the DC (ie, 
f=) term of the FFT to zero. The effect is nothing more 
than removing a DC level back in the signal domain, obdut 
Mellin transformation into the IM domain leaves the 
spectrum dependent upon the scale factor k [22]. Setting 
the f=0 coefficient to zero corresponds to setting h(f) to 


zero for ® f<1, where the unity upper limit is chosen 
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without loss of generality. The resulting transform of a 
scaled signal domain h(f/k) is 


wy. (/RIE AL = he? JR trae (35) 


which is obviously dependent upon k. The closer k is to 
unity, the smaller the effect. By increasing the f spectral 
resolution, the error can be reduced. The error may be 
insignificant for many applications (S-8], dut the 


technique should be used with care. 


The first solution has highlighted the need for a zero 
point correction. Another common solution [18-11] is 
developed by breaking the integral uo as before. Avain 
using unity as the upper limit of the left integral while 
remaining general in application, 


(hee eat a CAE eee 
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Two assumptions are made to get the correction term. First, 
that h(f) remains a constant h(@) over the interval of f 


between zero and one. Second, and not as easily accepted, 
ay Lp eeene 
oA 


‘oo 


= eA - O 


ae (37) 


Equation (37) pretty well shows why this assumption is 


suspect, but playing along for the momment, the question is 
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reserved ror a aver detailed look. Accepting the 


Seeurouron, the correction factor becomes, 





Lo) £00) 
Zes) = — a a (38) 


Because H(s) is a complex function, Z(s) must be applied 
(added to the imaginary part of the succeeding Fourier 
transform) before the magnitude is take to remove scaling 
@ependence. This correction is specifically derived for use 
with a continuous Fourier transform such as the optical 
Fourier in imaging systems with the added stipulation thet 
h(f) be nearly constant over the range @f<k where k is the 
largest scale factor expected. If an FFT is employed to 
rake this final transform, another correction should be 
applied as shown by Zwicke and Kiss [11] below. This 
@emrection factor differs from the first due to an 
invariant property of the FFT. The FFT of two unit step 
functions that vary only in scale are balanced. 


loa’ : * ; 
| Zc) 2 >. = ee |Z 0 . ¥ en? re wi e| Be 
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where mand p are arbitrary integers greater the zero and 
less that M. Successive FFT coefficients are summed, anéd 
the average value of the contributing terms taken resulting 


ine 


C= cot Crk /m) 


og 


fioaseis then multiplied by h(®) to arrive at the FFT Mellin 


eorrmection factor. 


A Co) : 
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When k/M is small the imaginary term dominates and the 
correction factor approaches that used in the continuous 
case (38). Most of the work done for the thesis on the 
exponential algorithm was done using the inappropriate zero 
point correction (38). Since its discovery was coincident 
with that of more powerful methods discussed in Chapter 


III, little data was taken using (41). 


To bound the error involed, an acceptable h(t) is 
defined, windowed and transformed using the Welling 
integral. The window limit is then allowed to grow 
unbounded and the resulting expressions are interpreted as 
the error. It’s deen assumed that for @<t<i, h(t)=h(d). 
Since the error resulting from the assumption in ecuation 
(37) arises from this interval alone, t>1 is ignored for 
the time being. Warping the signal as before h(x)=h(@) for 
all x<@ and h(x)=@ for all t>@. Evaluating the integral to 


a finite window width (T), 


or 
ae) Zoo} aia 
4 (oe) Jie wile ae) == (i-e ) 


Sin, CW my -30 Y% 
SS ee” Ge 
Cw 7/2 





=[r&col 


(42) 


4¢@ 





The term in brackets is the magnitude of the contribution 
from the region of integration, -T<x<d. Note that it 
involves a sin()/() term. The effect is to add a peak of 
(T)h(@) at the origin. The size of the peak depends 
directly on the window border T. Letting T approach 
infinity raises the spike at w =@, with lesser peaks at 
W=(2i+1)t%/T, where i is an integer. Fach of the subpeaks 
has a magnitude of 2h(0)T/((21+1)). Sudstituting W into 
the second relation yields the envelope (in brackets) and 


Phase as T tends to infinity. 


Bek (ome Se, 
Lev eee a [ Lo | (43, 


T@ 
The error bound in brackets, does not depend on the sample 
rate, or the size of the window. Any approximation 
approaching zero will have the same bound. Although the 
magnitude of the error is in a convenient form, the phase 
is indeterminate. For the correction to be applied, the 
complex addition must occur prior to the modulus being 
taken. This cannot be done, leaving the error uncorrected, 
but is bounded by 2h(0)/H for continuous and aperiodic 
discrete Fourier transforms. For the FFT, equation (43) 
does not bound the error. The sum of i/i does not converge 
as i tends to infinity. At any point on the FFT this sum is 
present due to the apparent folding. The error itself is 


not unbounded because phase differences in the sum of the 
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errors at any point may result in the envelopes adding 
destructively, reducing the actual error. Adding the 
envelopes is an unrealistic, worse casé approach. The FFT 
correction (41), can not be compared to (43) for these 
reasons. The other error correction, using the assumption 
in equation (37), can be compared. The first, setting h(@) 
to zero, or just plain ignoring the o<t<i interval have the 
error function bound in equation (43). Although only 
differring by a factor of two in magnitude, the constant 
phase is arbitrary, and equivalent to setting T=0;5 that is, 
assuming h(@)=0 over @<t<1l. This was the very problem the 
correction was developed to remedy, but is without effect. 
Equation (41), the correction for the FFT is not completely 
accepted by this author, albeit no real empirical evidence 
have served to verify or dispute the claim. Suspicions are 
raised on two aspects. The indeterminate phase of the error 
(43) in the continuous case arises naturally from 
approaching the t=@ limit with the Mellin intexvral. This 
quality is conspicuous by its absence in the FFT error 
correction. The FFT error correction is computed by summing 
the FFT coefficients in the complex plane. The average 
position of the resulting polynomial is the- error 
correction term. Eowever, the FFT coefficients of a finite 
duration signal are the values of the signal, evaluated at 


M evenly spaced points about the unit circle [3,23]. If 





the sequence is a constant, the average value is the origin 


ae une corplex plane. 


The zero point corrections for the Mellin transforr are 
unbounded at ™=@. More time could have been spent 
determining the best applied correction for the specific 
case at hand, but direct methods are developed in the next 
Chapter that obviate the need to employ the correction at 
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C. TESTS AND RESULTS 


It 415 worth admitting at this point that the results 
using the exponentially warped algorithm to achieve a 
discrete Mellin transform have not been good. More recently 
developed techniques in Chapter III greatly surpass the 
results reported in this s@ction. Although much of the 
theory used to improve performance in the following 
chapters could have been used here, this was a preliminary 
attempt that was later abandoned. The FM PEOces scl 
described in the previous section was built using FCRTRAN. 
Appendix Ais the documented program. This section wiil 
review the processing with actual plots of the signatures 
at different stages, discuss the required tests, introduce 
the testing approach, and finally intepret the results. The 


functional block diagram of the preprocessor, figure 5, is 
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near the beginning of this chapter and could be profitably 
reviewed. Figures 6 through 9 represent the step by step 
view of the signal processing, where the signal, a test 
Shape, is shown in figure ¢. The waveform, appears as an 
envelope, and is drawn with vertical lines indicating the 
sampled series. Figure 7 is a picture of the FFT, in its 
sampled version showing its symmetries. Figure @ is a very 
close approximation to the continuous periodic transform 
achieved by filling zeros to obtain the requisite spectral 
resolution. A + on the transform plot indicates an 
exponential sample point interpvolated from the sixteen 
unique points in figure 7. The warped samples are sent 
through the FFT block once more with the result shown as 
figure 9. Heavy spectral coloring by the h(2)/W correction 
factor is evident. Only the first half of the spectrur, 


from zero to the folding frequency, is valid. 


Complete scale invariance was never realized although 
its effect was greatly reduced. All the testing done on the 
exponential algorithm was an attempt at achieving and 
verifying shift and scale invariance. Much effort was spent 
structuring the tests to avoid the effects of processing 
noise so the actual algorithmic characteristics could be 
determined. Along the way, requirements arose to select a 
Suitable interpolation polynomial, an optimal zero point 


correction and other system improvements. In all cases, the 
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test procedure was to first select canonic test shapes. 
Squares and triangles were most frequently used, being 
shifted, scaled and combined to determine preprocessor 
characteristics. Scaling was most frequently by a factor of 
two or less. This corresponds to an aspect angle change of 
6@ degrees from the unscaled case. For many tests, care was 
B@ecatically taken to exactly scale a sampled series 
instead of the waveform envelope. The variation in input 
Signal when this isn’t done is evident when considering 
figure 10. For the envelope shown, the sampled series 
Cannot reconstruct the same signal, and may result in 
feature space variation. Two feature qualities were 
monitored in each test to determine preprocessor 
performance; insensitivity to shifting and scaling, and the 
ability to differentiate between canonic classes. These 
qualities were measured by visual comparison of the IM 
features, by computing the correlation coeffiecients and 
the mean squared error between the feature sets Di se 
differently scaled similar test shapes, and by computing 
the distribution of the error over the feature space. The 
latter test was an attempt to locate feature regions of 
class commonality, and regions of distinction between 
canonic classes. This approach allowed macroscopic and 
microscopic examination under changes of scale, shift and 


Shape. The clustering and separation qualities are data 
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dependent. Since no ship radar video data was used, all 
observations were with respect to canonic classes. No 


strong groupings were detected. 


As alluded to earlier, the results for the exponential 
algorithm were less than satisfactory. Test Shapes were 
carefully designed to minimize sampling effects, and to 
ensure low side lobes in the frequency domain. Many tests 
were conducted using each of the zero point correction 
methods with varied shapes, sample rates, and spectral 
resolutions. Classification on the besis of signal shape 
was very poor. The strongest correlation was between shapes 
of common duration. Table 1 shows a typical result of 
comparing a rectangular shape and a raised ramp. The 
scaling in each case was by 2 (€@6 degrees). Equation (38), 
Z=h(@)/W was used, but the others offerred Li ttice 
improvement. Consistently, the strongest similarity was 
Shown between shapes that had the same sample length, vice 


shape. 
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TABLE 1 


Canonic Shape Fourier - Mellin 


Feature Comparisons 


a. Peak Correlation Values 
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IIT. DIRECT MELLIN TRANSFORMS 


The last chapter develoved a method of obtaining the 
Mellin transform by exponentially warping the signal prior 
to using an FFT block. This technique is referred to as the 
fast Mellin transform (FMT). Although the promise of scale 
invariant features is attractive, some of the problems 
encountered that make the EMT unattractive are reviewed 
here. The required sample rate varies with respect to the 
data, making general applications difficult. The tendency 
is to use more samples than required, which quickly becomes 
costly in an exponential sampling schere. The need to 
exponentially warpv (to interpolate) a set of new samples, 
is expensive in time required. For true scale invariance, a 
correction factor is required but because of the integral’s 
unbounded nature at zero this correction factor is 
indeterminant. several correction Methods have been 
employed, but they depend on unspecified date 
characteristics. These effects combine to make the actual 
performance of the algorithm poor. Although scaling erfects 
are mitigated, they remain an artifact, which is disturbing 
to classification attempts. This chapter outlines the 
effort to remove these limitations. Some useful Mellin 


properties are developed, and then applied to establish 
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several Direct Mellin Transforms (IMTs) which were built, 


and their performance compared. 


A. SOME USEFUL PROPERTIES 


Some general observations are made here about the Mellin 
transforms, and are followed by some specific relationships 
which were derived and applied. A property of the Mellin 
Ss-domain is that it is unaffected by scaling changes in the 
original x-domain. Figure ila is a test shape in the 
Xx-domain. Two features are identified according to their 
amplitudes A and B, at x=a and x=bd respectively. The ratio 
of a/b equals c. To be simply scaled by k, h(x/k) must 
remain the same in all aspects except that the distance 
between features has been changed according to the scale 
factor k. Figure 11b shows the scaled domain. Features A 
and B are again identified at their scaled positions ka and 
kb. The signal property maintained by scaling is relative 
positional integrity, that is, ka/kb equals c as before. 
The positional integrity of the features, the ratio of 
their distance from the origin, to that of another’s is 
unchanged. Restated, an operation Ofh(x)] in the x-domain, 
will leave the s-domain modulus invariant to simple scaling 
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Mabe tnat the entire domain is scaled, so no unscaled shift 
in the domain can be permitted as already discussed. In the 
Mellin s-dorain, any simple scaling in the x-domain results 
in a phase distortion (the modulus is invariant to k). 
Manipulations in the s-domain will leave that domain 
invariant to k, as long as the modulus is modified by a 
multiplicative factor of constant phase. That is, if {G(s)! 
is an arbitrary function of s (except that it does not 
depend on k) and {M{h(x/k)] /=/H(s)] is also invariant to k, 
then their product is invariant to k and the x-domain 
remains simply scaled. For instance, in the x-domain the 


operator Ofh(x)J=x(h(x)), does not meet condition (44). 
ab(apfha) F FCxZ&) (45) 


So the Mellin transform’s modulus of xh(x/k) cannot de 


invariant to k. 


In Chapter II, a error was apparent because h(x) was not 
equal to zero for x=@. If h(x) could be modified with an 


operation that met the condition of equation (44) and 
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always produced a series that was zero at x=8 a general 


approach can be developed. Consider two operators, 


Of hx7e) ] = F(x AC x/a)) £ (x/h) 


OL &¢«/a)] = « 4 (Acx/&)) = £¢x/h) a 


Equation (4€) will produce an acceptable f(x) as df/dx=@ at 
x=@. Equation (47) will always produce the required 
condition. To see the frequency domain equivalents, we must 
S@eoume that the Mellin integral exists. Integrating by 


parts, 


J dewra) eee ~ | Hes) | 
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The result in the frequency domain is consistant with the 
conditions stated above. The limit of h(x) as x goes to 


zero or to infinity must be zero if equation (48) is to be 
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true. Similarly, under the same conditions, it can be shown 


that 


[M[ cs hexihy) | [ (se Hes | (49) 


Equations (48) and (49) clearly do not apply where h(@) is 
not equal to zero. However, by using the x(ad(h(x))/dx) 
modifying operation of (47), a series can always be zero at 
x=@. The function h(x) is further constrained by the fact 
that it must fall off faster than 1/x. This assumption must 
be valid and the modifying operator must be applied in the 
x-domain for the Mellin integral to exist in general. A 
Mellin transform of a function, after having a modifier 
applied, will be called a modified Mellin transform. 


= © ot £ex) ig 
Ha ¢s) = \ Ve 1X Ax (58) 


The integral is close to a form which is realizable, except 
for the upper limit. For a finite sampled series, h(n) will 
be assumed zero outside of the interval @2n<N. This 
truncation effects the transform of an otherwise infinite 
series. In this application the Mellin is applied to an FIT 
frequency spectrum. The truncation in the frequency domain 
is due to band limiting the signal prior to sampling. 
Scaling in the time domain will not result in simple 
scaling, but will adda dependence on the scale factor xk 


that can not be removed by the transform. An approach by 
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Prost and Goutte is used to predict the size of the error 
[24,25]. First a suitable function will be selected and a 
relative error of truncation (RET) determined and applied 
to two scalings. The relative difference of the feature 
wees is found and identified as the error. Remerbering 
that this is applied to a frequency spectrum, dh(f)/df is 
Seproxri@mated by a function of the form. 


AS GR =n ad 


ane (o2) 


The modified Mellin transform of (51) over a finite 
frequency range would te approximately, 
Ha *¢s,F) = mai mctMe hil 

Ye (S2) 
The lower limit has been set in a manner to be consistant 
with Plancherel’s theorem. A convenient worse case 
assumption is that the lower limit is essentially zero, but 
this depends, in general, on F and the data itself. 
Ha*(s,f) converges toward Ha(s), equation (£8), in the mean 
square as F tends to infinity. The mean square error will 
not provide an expression for the error that can be used to 
correct for it but is useful to measure the effect of the 
truncation in more general terms. A relative error of 
truncation (RET) can be computed if the error is assumed to 


be distributed evenly over the range from which Ha*(s,f) is 
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computed. By employing Plancherel’s theorem, the mean 


Squared values are, 


eC, = iceta ae jo fr er** at ae.) 
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(54 ) 
The RET is defined and solved as, 
TS ue a 
RET = (a | = € (2F*+F +1)? (55) 


But the limit F depends not only on the pass band F, but oa 
the scaling in the f-domain. A relative error (RE) between 
a truncated spectrum and a scaled and truncated signal 


would be more complex, but worth the effort. 


e“(2 F* +F +!) -( l/h} Paarl 2 bear To l) 


RE = 
[newer et) 


(56) 
where @<k<1. Two observations should be made here. First, 
the relative error of truncation (5&8) and the relative 
error or difference between two truncated spectrums scaled 
differently (56) both depend on F. F depends on the cut off 
frequency of the low pass filter prior to any sampling. It 


is often chosen to minimize aliasing depending upon the 
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[Merete tions of the Sampling circuit. Second, RE depends on 
k as well. The importance of scaling differences to this 
data type can be readily seen. If equations (55) and (56) 
are valid, and if the range of k can be bounded (a design 
specification) F can be chosen to realize a stated RE. Or 
if F is fixed and k bounded, the RE may be determined for 
evaluation. If the data type is not apvropriate, a more 
representative function may be determined and used in place 
of equation (51) to attain better expressions for FET and 


RE. 


3. ALGORITHM DEVELOPMENT 


Part A above provided some background for making the 
Tirect Mellin Transforms (DMTs) in this section. In this 
presentation the justification is given with the 
application forest. However, chronologically the neat 
package presented above was preceded by extensive 
evaluation of empirical results, not vise versa. Appendiz 3 
documents the FORTRAN implementation of all the algorithms 
developed in parts 1 and 2 below. Figure i2 is a functional 
block diagram of a generalized FM preprocessor using a DMT 
to show the simplicity with which it can be applied, as 


opposed to the FMT covered in the second chapter, figure 5. 
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1. First Difference Approximations 
Although developed through a different rationale, 
this first algorithm was developed by Zwicke and Kiss [11]. 
starting with a sampled series hi, the series is operated 
on by the modifier defined in equation (51), using the 
first backward difference to approximate the derivative 


with respect to x. Unit step size is assumed. 


x A Kx) = m(k,-d _jamVA, (eis, 
d x ma 


Taking the trapezoidal rule to evaluate the modified Mellin 
integral (52) while recalling that h(@)=0, and h(N) is 


assumed zero, 


N-/ 
Ha, (™~) = > vA, mm? 


(58) 
na 
where s=-j2mm/M. The complex ccefficients are 
miz cos 20 (am/m) Dam sina (om Sa) La im) lee) 


They can be calculated off line and stored to produce just 
the desired characteristics. The factor (2#m/M) could be 
any number that produces an interesting feature. Undesired 
features need not be computed (what in general was aN by M 
Watrix, where ™ is the number of Mellin transform 
coefficients and N the number of x sample points). If a 


relatively small number of features is required, perhaps 
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the processing will be manageable. Notice that no zero 
Meant correction is required. Cnly data changes contribute 
to the transform. These observations are valid for any of 
the modified Direct Mellin Transforms developed in this 


section. 


By using a central difference instead of the 


backward difference, a similiar result is obtained. 


N-! 
Har Com) = » CA, ~-4,_, ) n° /2 (69) 
ms | 


Other numerical integrations May be used with improved 
results, and other methods can be used to increase the 


order of the approximation. 


To test the algorithms, a ramp and inverse ramp were 
used in a scaled and unscaled mode. The ramps and their 
Scalings are shown in figure 13. Figure 14 is the analytic 
results of both waveforms plotted with the transform found 
using equation (€@). This is a dramatic improvement over 
anything used with the methods discussed in the previous 
chapter. The noise of the signal appears to be diverging as 
frequency grows, but over Lhemnene plotted, the appearance 


is that of an algorithm trying to do well. 
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4 (2) 


b§ TRIANGLE 


Test Shapes 
Figure 15 
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<. Second Llifference Approximations 


A second difference algorithm can be achieved by 
proceeding as before. The pure Scaling operator used _ to 


prepare the signal is, 


G 
0 ike a 2 of Acx) — ~ + 
C a) | ve = neh, DIN OL Ree 
And the new algorithm is, 
; 
. - st 
Hew “5? > 5iay D. ath, (E2) 


me | 
Other second order operators nave been used, and 
partitioned into forward and backward difference variations 
to (61), but this appears to be a basic and useful form. 
The color i/(sti) is present to approximate the modified 
Mellin of equation (68). The term (sti) is valid assuming 
that dh(x)/dx is exclusively upper bounded by In(x)/x as it 
approaches zero or infinity. For comparison, another second 
difference algorithm was developed based on the modifier 
acwex )i(x(d/dx)h(x). 
M-1 
Hoc DS math eld, - £1) ow 


S 


may 


(63) 


morse i1s rowehly the sur of the methods defined in equations 
(68) and (62) above. The assumption for deriving the color 
1/s is even less restrictive than before, but the algorithm 


performs poorly compared to (62). These transforms depend 
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eseGonc Oirference characteristics, but are modified by a 
1/s term which has a stabilizing effect. This is equivalent 
momandivision by Xx and integration in the x-domain. The 
order of the approximation has been increased by the 
modifier. Results using equation (68) and (62) should bde 
alike. Figure 15 is the results of using (62) compared to 
mmemenosedc form solution to the figure 13 test shape 
transforms. An improved performance over (€0) is seen over 
some of the range, but a drop as frequency increases 


degrades the accuracy in figure 15b. 


5. Higher LTifference Approximations 


Higher difference approximations can be developed. 
For instance, one algorithm depending upon the third 


difference is, 


- wiligs 3 s*2% 
Ha; = sé > vA. m haa. 


The performance of higher order algorithms becomes 
increasingly suspect because of the extreme weighting they 
apply to different parts of the series. Different 
algorithms exist, but this weighting is always a factor. A 
smoother transform is achieved, but a large érror is likely 
to develop due to the algorithm’s dependence on higher 


order derivatives and the nature orf sampled data. Eowever, 
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these comments are speculative since they were not 


confirmed experimentally. 


4. Higher Order Integrations 

Higher order integration rules should be able to te 
used with a corresponding improvement in performance. One 
using the first difference with Simpson’s rule Was 
implemented with dissappointing results. Figure i6 is the 
result of such an implementation. The droop for the ramp 
input is apparent even though not present in the 
trapezoidal rule used in subsection 1 above. The higher 
frequency error is also more prevalent than before. A 


pmeeram E€rror is of course suspected, but was never found. 
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IV. CLASSIFICATION PREPROCESSING 


The problem of determining important Signal 
characterisics can be approached in at least two different 
ways. First, by trying to learn what is important in human 
me@orni tion, and then trying to adapt a machine to emulate 
that behavior. Or second, by using successive transforms to 
rerove inno rina t) On known to be superfluous to 
classification, while keeping enough it 0 bina tO to 
reliably assign an object to a class. Addressing the 
hommenr, even though it is difficult to determine specific 
details, some key aspects of human visual recognition are 
discernable. Chief among these is that the intensity level 
of a scene, or object, does not appear to be as important 
as the relative position of the edges, i.e. the shape 
separating different intensities and frequencies [25-28]. 
Examples in scene analysis show clearly that the edges or 
Shapes are far more critical to human recognition than the 
relative power differences themselves. The invariant shapes 
or angles in scene analysis find their analog in ratios 
between similiar points in differently scaled time series. 
Examination of many range only radar video ship signatures 
has provided a basis for noting that the relative position 
of time domain features remain constant over changes in 


aspect angle, while the relative intensity of the features 


Ve 
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vary greatly as shown in figure 17. As the aspect changes, 
Mmemvars Set Of ratios remains constant within a ship class, 
then the set is identified as information worth preserving 
for the classifier. Conversely, since relative intensitiy 
is not a stable measure over aspéct angle, or from among 
different ships of one class, that information should be 
intentionally removed to provide tighter natural 


clustering, with the minimum number of features. 


A. INFORMATION REQUIRED TO CLASSIFY 


There are two preconditions that must both exist for a 
set of possible input signatures to separate into distinct 
classes. First, the features of a particular class must 
have some common characteristic about them, and second, 
this characteristic must in some way be unique with respect 
to other classes. The assumption in existing radar 
Signature classification projects is that there is enough 
information Tei the Signatures to permit this 
Slacoitfication. Short of actWalyetrying to classify with a 
Seve Ot redlistic signatures, stue anmelytical determination 
that sufficient information is present in a set of all 


possible ship signatures is difficult to approach. 


To establish how well the autonomous classifier is 


performing some measure of the classifiability of the set 
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received signals is desired. Failing this, some discussion 
of the information capacity of the preprocessor should at 
least be considered. The preprocessor produces the FFT 
Tragnitude of a sampled signal as the output of the first 
stage. Most of the unique positional relationships of the 
Signal upon which human recognition apparently depend has 
been destroyed. Next, the Mellin transfrom stage 
effectively distorts the signal and uses the magnitude of a 
second Fourier transform as the output features. Signals 
reconstructed on the basis of FFT phase information alone 
usually provide sufficient similarity to be associated with 
the original signal, whereas reconstruction on the basis of 
Tagnitude does not retain any Significant detail except 
when the signal is symmetric [29]. The data is also masked. 
For most applications, the data is frequency shifted, 
filtered and sampled as a baseband signal. This windowing 
Tasks the magnitude characteristic and is specifically 
designed into the processing. The resulting FM features are 
insensitive to positional and scaling relationships that 


are necessary in human visual recognition. 


allvt vc Support is also available to aquantify the 
importance of relative position of events [29]. By 
considering rms error (due to spectral phase and amplitude 
quantization for random signals) it has been concluded that 


approximately two more bits are required Orr the 
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quantization of phase information than amplitude for the 
same rms error. A separate analysis applied distortion rate 
mmeory tO real~part, iraginary- part, and magnitude-phase 
encoding of the DFT of random sequences. The result was 
that phase required 1.4 bits more storage than magnitude 
for a similar error [30]. A third approach concluded that 
the Fourier phase includes 1.8 bits more information than 
the magnitude [31]. This was based on analysis of image 
reconstruction from kinoforms (phase-only holograms). The 
fact that phase-only reconstruction preserves much of the 
correlation between signals would suggest that the location 
of events tends to de preserved. Further, it seems’ that 
this information is lost by taking only the magnitude of 
the Fourier bean st OR Ms Another Pnveres tite. albeit 
informal, view is apparent as one considers the phase-only 


Signal as a spectral whitening process. 


FL fe] = F(f) = | Fc#)| ang (65) 


For reconstruction by phase alone, where the magnitude is 


set to one; 


(66) 


FL 4, 41] = aa FL fee] 


Since the received radar signature will have an abundance 
of low frequency spectral lines and smaller high frequency 


components, the low frequency information is not weighted 
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as heavily as the higher frequency information in the phase 
reconstructed signal. This seems like it would accentuate 
sharp changes in the reconstructed signature without 
removing relative location information. The result is the 
summation of the different frequency components, all with 
zero phase (i.e., no positional or amplitude distribution 


information can remain). 


Considering the Mellin transform next, in a continuous 
case, the exponential warp does not lose any information. 
To zero the first data sample as required by the Chapter 
III modifications, will surely destroy information, but 
this may be confined to the [TC term alone. The information 
lost during the final transform and magnitude is difficult 
to assess. Increased masking occurs due to the spectral 
truncation, so actual information loss may not be as great, 
but masking distortion may be greater than before. So 
approximately two bits of information are lost. Only a 
quarter of what was, remains. Information is also lost when 
the transforms are normalized in the processing so that any 
power calculation is also meaningless. Some interesting 
questions arise. After removing positional relationships, 
Scaling, and power, what signal qualities remain and are 
they useful in classification? Although it is true that 


this insensitivity may add a certain robustness to the 
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syster, the arbitrary loscemort Valid “classification 


ieeor@mation showld be rinimized. 


By exarining the quality that must be ignored, and bdy 
comparing its removal to what is actually removed by the 
processing, an interesting result will develop. The effect 
Sree tire domain shift on the frequency domain is an 
additive phase term, linearly related to the frequncy of 
the coefficient, as in equation (2). Most of the structure 
of the signal is held in the phase relationships with 
respect to the fundamental and higher frequency terms. So 
shift can be defined as the phase of the fundamental 
complex coefficient. sy setting the phase to zero, and 
adjusting the other coefficients according to their 
GComponemt frequency, the structure of the signal is not 
lost but reconstructed about the fundamental as before. If 
there are N/2 spectral phase angles, only the fundamental 
needs to be zeroed to remove the shift. If the information 
is contained uniformly in the spectral phase, then only 2/N 
of this structural information needs to be removed. When 
the magnitude is taken to produce shift invarient features, 
all of the phase relationships are destroyed. {(N-2)/N of 
the information once held in the phase was removed 
needlessly. The amount of information lost removing the 
shift can be made arbitrarily small. As N grows unbounded, 


the amount of information that needs to be removed tends to 
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zero. This result seems to be supported by experience. With 
enough samples of a shape, its position with respect to the 
observation field is immaterial. In theory the principle is 
sound, but some practical limitations may degrade predicted 
performance. Recalling that an exponential warp translates 


scaling to shifting generalizes the result a bit further. 


Ait = # (er A) (67) 


The same priciple that permits simple shift removal is also 
valid for the removal of scaling dependence as well. Using 
the Mellin transform, scaling dependence may be removed by 
zeroing the fundamental and adjusting all the other 
coefficients as described above. For the FM preprocessor, 
the information lost removing the scaling and shifting 
dependence may be made arbitrarily small by increasing the 
number of spectral samples used. Since the number of 
spectral samples can be increased by filling zeros onto the 
finite signal in the original domain, this process does not 
effect the data sample rate. This approach was not verified 
experimentally, but represents a potentially powerful tool 


to analyse and improve the extracted feature space. 
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B. RANGE CNLY RADAR 


This se€ction addresses ship classification on the basis 
of information gathered from a range only radar video ship 
Signature. An example of such a signature has already been 
considered as figure 16. Classification by range only radar 
Signatures is subject to the same distortions discussed 
above. Typically, the radar return is detected and isolated 
in ae range gate that is sampled and digitized. The 
rectangular sampling window can be considered the range 
gate itself. The range gate is designed to ensure that the 
included range is greater than the maximum ship length so 
that the time/range windowing has no effect on the 
frequency spectrum of the signature, other than increased 
spectral resolution. The placement of the ship signature in 
the window is not set, neither from encounter to encounter, 
nor from pulse to pulse (jitter). The total effect joins 
together to produce the framing distortions. Signature 
Scaling results from viewing the ship from different aspect 
angles. The sampling rate must be done at more that twice 
the imverse of the resolution of the receiver. Quantization 
levels are chosen in a manner to reduce that predictable 
random noise to an acceptable level. The pulse to pulse 
jitter is an ever present characteristic of the radar 
problem, but at a normal resolution (greater than 25& feet) 


integrating the return in the time domain removes’ the 
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[eeemem owen tect, reduces scintillation, and improves the 
resolution of the Signature. Although predetection 
(coherent) integration is more efficient, post detection 
(noncoherent) integration is more common because of the 
convenience of not having to preserve the radar frequency 
(RF) phase. For post detection integration of n pulses, the 
Signal to noise ratio would be something less than n_ times 
the signal to noise ratio for one pulse (32]. More 
important to the recognition problem itself, for a stable 
system by the law of large numbers [33], fluctuation of the 
average value of the return will be overcome. That is, with 
the integration of n pulses the resolution (R) will become 


finer as 


R = | S* Dm (68 ) 


with a probability of 1-D, where S squared is the variance 
of the Signal from pulse to pulse. For very high 
resolution, the cost of making the signature stable with 
respect to the integrator oecomes prohibitive. It has 
become convenient to integrate the spectrum of the 
Signature bécause the jitter effects can be completely 
removed. Another limiting factor in integrating over a 
period of time is that the position and aspect cf the 
target are dynamic. They change with time. A conceptually 


attractive solution is a recursive filter which weights the 


81 





integrated pulses such that the older they become, the less 


weight is accorded to them. 


If the course and speed of the target ship are known 
from measurement of the target track, it is possible to 
infer the aspect of the ship. The range profile can give 
some estimation of size. Unfortunately, the three 
dimensional change in aspect angle, commonly suffered by a 
Ship presents more than just a video signature scaling 
change. The radar cross section of even individual 
Structural components of the radar target changes with 
respect to the aspect angle. The composite effect is that 
Ship signatures vary greatly with aspect angle. The radar 
is an electromagnetic sensor, reacting to energy reflected 
from the target. These reflections are a result Of 
Scatterers that are related in dimension to the wavelength 
of the illuminating energy. Because of the great difference 
in wavelength between light and microwaves, what can be 
"seen by radar may be quite different than that seen by an 
eye. Also, when measuring size or any distances with a 
radar of high resolution (less than £@ feet), an error can 
be incurred since the extremities of the target are not 
always good scatterers. Echoes from the forward or stern 
portions of the target might be observed in the noise, 
especially for a relatively low power radar [32]. After 


reviewing hundreds of signatures, it appears that major 
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features, such as overall ship length and dominant mast 
Suructure are fmequently discernable, but vary in relative 
amplitude. Resonance, shadowing (one reflector hiding 
another), multipath returns, the mapping of three 
dimensional aspect changes onto a one dimensional time 
series, and the amplitude and phase of component returns 
Summing constructively or destructively to cause a 
Scintillation of the composite target. Some of the 
variations caused by these conditions can be lessened by 
integration but major effects remain causing the signature 
to vary in shape and content with the aspect angle. for 
this reason, a class feature volume cannot be reduced to a 
Single point, but will remain a hypervolume in the feature 
Space even in the ideal case. Any selection of features 
Should try to minimize this volume. Features should be 
selected that are relatively insensitive to known 


superfluous effects. 


C. CLASS DISCRIMINATION 


In the last chapter, the major concern was removing two 
sources of variance with no classifying value. Algorithms 
were developed and canonic shapes generated to verify the 
algorithms and demonstrate the invariance to scaling. In 
the same manner, this chapter has reviewed information loss 


and process masking. The question of whether sufficient 
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PorOrmation is present to classify was raisec. To answer 
this question some simple classes of canonic figures are 
defined, and put through the entire FM preprocessor 
documented in Appendix FE. The DMT algorithm found to offer 
the best scale factor rejection in Chapter III was used to 
generate the final features. The algorithm chosen is based 
Oma second difference modification and is defined in 


equation (€2), 


N-/ 
| 
Haz (3S) = Say >, Arh ont?! (62) 
m= | 


After canonic tests were made, preprocessor performance on 
several ship signatures was recorded. Although this was 
premature in the logical test sequence, the results are of 


some interest. 


1. Test Shapes and Results 

Four test shapes were used. Figure 18 shows the test 
Shapes. All were scaled and shifted originally to test for 
algorithm verification and demonstrate scale invariance. In 
this series of tests they were left fixed and used in 
different combinations to try to detect shape presence in 
the FM feature space. The object is to differentiate 
between different canonic classes. Figure 19 shows a 
comparison of the shapes, rectangle and triangle. A _ test 


combining the rectangle and triangle was the subject of 
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figure 29. Finally in figure 41 a rectangle with two 
triangles is shown in the IM feature space. It is clear 
from the plots that most of the scale variance has been 
removed. Just as important, some quality does remain that 
differentiates between the canonic shapes. A Square ia 
general can de differentiated from a triangle. A ship 
with a single mass of superstructure can be separated from 


one with two such masses. 


Although it’s clear from the plots that there is a 
unique quality left in the feature space to allow the time 
domain shapes to be classified, some quantified measure of 
System performance is required. The magnitude of the 
feeature vectors have all been normalized with respect to 
Geet irst coefficient, sO 1L)a that region little 
discrimination can be expected. For higher Mellin 
frequencies, noise dominates. A region of coefficients, 
11-1908 was chosen to classify the shapes by correlation. 
The results are included as Table 2. The improvement in 
performance over that shown by Table 1 in Chapter II is 
dramatic. The methods in that earlier test resulted in the 
observed length of an object being the distinguishing 
Speeceni1on for classification. Using methods supported ia 
Chapter III, variations due to scaling and shifting of the 
original domain have been removed. The features now reflect 


the shape of the object in the time domain. An unusual 
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effect is noted with respect to the difference squared 
analysis in Table 2b. Although TRI2 is more closely 
correlated to TRI1 than RECT2 the squared error shows just 
the reverse. The results are encouraging. The preprocessor 
has greatly simplified the classification problem for the 


canonic shapes above. 


e. ohip Signatures and Results 


A single ship was used to make these preliminary 
tests for this thesis. Signatures were taken every ten 
degrees around a ship from zero to fifty degrees. The 
results are plotted and compared over twenty degree aspects 
in figures 22-23. The signatures are the result of very 
high resolution radar signature data that has been degraded 
and smoothed to a lower resolution with essentially no 
noise present. Recalling that the purpose of the 
preprocessor was to remove variance due to pure shifting 
and pure scaling, leaving enough ini@ortration for 
classification, to the eye there seems to de little 
encouragement from these results. It is recalled that a 
goal of this preprocessor is to make the classifiers job 
easier by removing dependence on shifting and scaling of 
the original data. The information that remains depends on 
unspecified signal characteristics that here appear to 


useful in discriminating shape classes and possibly ship 
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classes. However no real conclusion can be drawn at this 
point because of the small data base and the atsence of an 


automatic classifier to generate an optimal feature Space. 
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Table 2 


Canonic Shape Fourier - Mellin 


Feature Comparisons 


a. Peak Correlation Values 
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V. CONCLUSIONS 


Mieee preprocessor design began by considering a generic 
classification system. A distinction was drawn between the 
classifier and the preprocessor. The preprocessor is 
meomeem Specific. It assists in the classification by 
extracting a set of features for the classifier. The 
extracted feature space has enhanced natural clustering on 
the basis of shape by removing information that was 
extraneous for classification. Two useless cnaracteristics 
were identified as shifting and scaling. The preprocessor 
was designed to remove dependence on these two 
characteristics by using the invariant properties of a 
Fourier transform followed in series by a Mellin transform. 
The resulting set of coefficients are Fourier-Mellin (FM) 


features. 


A. REVIEW 


In Chapter II a Mellin transform was developed using the 
conventional Gu 2aa a )! processing approach which 
exponentially warps the domain and then transforms the 
Spectrum by an FFT. This method is sometimes identified as 
a fast Mellin transform (FMT). The unbounded behavior of 
the exponentially warped frequency spectrum was shown to 
fesult in a funetion that could not be transforfned. That 
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is, the warped function has no transform because the Mellin 
integral PSeeeaeeterrimate “at = the lower limit. ~Lrror 
correcting techniques cannot compensate for the effect 
because the error itself cannot be computed in general. The 
bound for the error was found and seen to be the envelope 


memeexisting error correction functions. 


Chapter III developed some useful properties of the 
Mellin transform that were used to modify the signal so 
that the pitfalls isolated in Chapter II were avoided. This 
was done by modifying the input to the Mellin transforr to 
always be transformable. To simplify the implementation and 
to control the effects of sampling, a direct Mellin 
transform was used for the developrent of the modifiers. 
Several suitable modifiers were determined and tested with 
differently scaled inputs for which the closed form 
Solution was known. The direct Mellin algorithm that 
produced features closest to the closed form solution was 
chosen for use in the preprocessor. With the modifications 
in place, the preprocessor was tested and shown to produce 
features that were invariant to shifting and scaling. It 
was also shown that the features retained enough 


information to classify canonic shapes. 


Chapter IV discussed what type of information is 


required for classification. Signal structure or shape was 
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identified as key information. A discussion of what 
formation is required for classification and a means of 
keeping the signal SeRUcCUCUreG intact throughout the 


preprocessor was advanced, but not empirically verified. 


B. FUTURE WORK 


The design IM preprocessor does produce a feature space 
with enhanced clustering, but problems remain to de 
resolved before the full potential of the system can ode 
realized. There are three extant conditions that detract 
from the performance of the implemented preprocessor. 
feos, the preprocessor 1s not computationally efficient. 
Second, most of the signal structure that should be vital 
to classification is obviously lost. And third, a complete 
verification of performance has not been conducted. The 
current preprocessor produces an enhanced feature base for 
a classifier, but attention to these main weak points will 


greatly impreve the applied techniques. 


eet Lol ent Brocessnme 


A direct Mellin transform using N spectral samples 
to transform into M Mellin spectral coefficients requires M 
by N multiplications. If only a few coefficients are to be 
used then the number of multiplications may be small. Using 


the FMT requires about N(2@+ln(N)) arithmetic operations 
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Pometne interpolation and the FFT. If twenty-five or more 
Mellin coefficients are required, the FMT is faster. The 
modifiers uséd to prepare the spectrum for the direct 
Mellin transform will also work for as FMT, but this should 
be demonstrated experimentally. Because the FMT directly 
weights the lower numbered samples it may produce more 


accurate results as well. 


2. Conservation of Information 

Chapter IV discussed the type of information 
required for visual pattern recognition in humans. The 
DPeeservation of this information should be a specified 
design goal for the preprocessor. The preprocessors built 
or this thesis removed much of these vital signal 
properies. A means was introduced to limit or control the 
loss of structural detail by zeroing the fundamental phase 
mmgedausting Gach of the remaining complex coefficients to 
meceemetruct phase relationships. This may be done in the 
frequency domain as described, or by a similar operation in 
the time domain, shifting the centroid to zero. In either 
case, to transform more information about the signal will 
efrectively increase the sensitivity of the features to 
Gi@ard@otenistics in the time domain. This increase in 
sensitivity needs measured to confirm the approach. It is 


€also possible that continued selective zeroing of 


a8 








coefficients may offer improved performance or robustness 


to the system as a whole. 


Seeeveriuticavion 

Although the improved performance was demonstrated 
with respect to ealier FM digital preprocessing, a direct 
improvement factor needs to be established. Time domain 
correlation should be used as a measure of original signal 
@eaoslriability. This weasure needs to be compared to the 
FM domain correlation recorded as Table 2 in Chapter IV. 
Next, the preprocessor or an improved version will have to 
be married to a classifier and realistic data used to 
evaluate its effect on the classification system. The 
preprocessor built was design to be used on-line. An 


on-line classifier needs to be built as well. 


The systematic evaluation of ship profiles using IM 
features is still a requirement. For several ships, IM 
features must be singled out and plotted together as a 
function of aspect angle. The purposes are to establish a 
range of aspect over which classification may be possible, 
to evaluate changes of structural content as discussed in 
section two above, and to determine the beam signature as a 
classification “node. Although these purposes assist in 
the system design, the third is more of an operational 


necessity. All ship signatures will degrade to the beam 


ee 





aspect node so that this hypervolume in the feature space 
is occupied in common. Therefore the bear condition rust be 
detected and withheld from ever entering the classifier. 
Nor should beam signatures be used for classifier training. 
The beam “node needs to be determined separate from the 
classifier. Initially this information can be plotted to 
examine feature behavior and confirm the rmrethodology, 
automated methods will quickly follow. These analysis 
functions may never reside in the classification syster 
itself, but must be a part of the tools used to develop a 


workable classification system. 


120 





AAA AAA NNAANANANANANAANANAANANANANA 


AAA 


150 
200 


APPENDIX A 


RKEVMEKKKAKAKKKUKASRESKKKKKKKKLRSCELAKEKAKKEKKKKKAKKEKKE 
* THIS FROGRAM IS a DIGITAL L[NPLENENTATION O * 
* A FAST FOURIER TRANSFORM FOLLOWED BY A MELLIN ~~ 
* TRANSFORM. THE DIGITAL IMPLEMENTATION Of THE * 
* MELLIN IS AS AN EXPONENTIALLY SAUPLED SPECTRUM * 
* THEN RUN THROUGH AN FFT AND CORRECTED FOR THE * 
* ERROR BY ANY ONE OF SEVERAL CORRECTION * 
* SUBROUTINES, CORCTX.. * 
VT TELL TTTST TTT TTT Tete TTT TT STE LT tel ete rere rere rer eS 
RECORDED PLOTS / PLOT NUMBER USING RECORD CALLS 

-TIME FUNCTION Neeerh yf 1 

~SAMPLED TIME FUNCTION / 1 

=Fre (CORT) / 2 

~EXFONENTIALLY SAMPLED FFT / 2 

~UHIFORMLY SAMPLED FFT Z 5 

~DISCRETE MELLIN ®EATURES / 4 


RECORDED PLOTS / PLOT NUMBER USING RECANG CALLS 
~TIME FUNCTION (CONT) / 1 
~FFT UNIFORMLY SAMPLED / 2 (MAG & PHASE) 
-MELLIN FEATURES / 3 (MAG & PHASE) 


N XBEAL (200) ,XI MAG (200) ,XINT (200) ,T (64), 
Sa 


U 
AMP (XREAL, XIMAG,N) 

ECORD (KREAL,XIMAG,N,1ISCALE) 
ET 

DY 


Se ie SAMPLze POINTS AS PROVIDED BY THE 


Hon 


u 
AMP (XEEAL XIMAG/ N) 
ECORD (XREAL, XIMAG,N,ISCALS) 


ow | 
tH 
Paes ok 
i) i | 


¢ OOH >= Il 


OWOSSO HMO RH ONMs NNO ONS Ste eANo 
Oe 


BH “HOM 


Git Il" 
rd il Obs Iie 
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310 


© 
© 


ANA = AAANAW 


© 
© 


ANANANANANNAS 


a9 


QO1 


ooh woh 
Ou 


ANANDA AnANSES 


Q 
O 


CALL FFT (XINT,PRT,L oN, 
CALL FFT XREAL, LIAAG, N, 


CALL RECOED ThA6 LESCALE) 
CALL RECOR XHEAL MEET "45 ne 

ISCALE = 3 

DO 300 I=1 | 

XREAL (I =S6aT (XREAL (I) **#24+X2IMAG (I) **2) 
XIMAG {I} =0. 

CONTINU 


CALCULATE THE NEW SAMPLE POINTS AND INTERPOLATE 
TC FIND THE EXPONENTIALLY SAMPLED VALUES. 


CALL NUPTS (XINT,M,N) 


AFTER NEW INTERVALS CALCULATED AND STORED TEMPORARILY 
IN VECTOR XINT, STORE FOR LATER PRINTLING LN 2 VECTOR. 
DO 400 I=1,M 

eat, 

a THe 

WITH THE NEW TIMES IN VECTOR XINT, AND THE CURRENT 
SPECTRUM SAMPLES IN VECTCR XREAL, COMPUTE THE NEW 
EXPONENTIAL SAMPLES AND ENTER THESE INTO VECTOR 
XINT BY USING THE CHOSEN INTERPOLATION METHOD. 
FINALLY, RECORD THE FIRST SAMPLE TO BE USED LATER 
TO CORRECT THE MELLIN TRANSFORM FOR LOW FREQUENCY 
LOST IN THIS EXPONENTIALLY SAMPLED TECHNIQUE. 
FO=XREAL 1) 

CONTINUE 

CALL lt all ladle 

WRITE (2 40 M 

FORMA (£4) 

CALL : AX XINT, PRT, M,SCALE) 

DO 410 

XREAL I}= =xtnr (I) 

XINT ( a SCALE XINT (I) 

XIMAG (I) =0 

WRITE (2 4 ASPET AT) ¢ XINT (I) 

FCRMA £10 F10.5) 

CONTIN 


SUBMIT THESE EXPONENTIALLY SAMPLED VALUES TO THE 
FINAL FFT BLOCK. 


CALL FFT (XREAL,XIMAG,M, 40) 


APPLY THE CORRECTION TERM, FIND THE MAGNITUDE 
OF THE NOW SCALE AND TIME INVARIENT FEATURES. 


CALL CORCT1 (XREAL, aN Pigcgl8 
CALL CORCT2 XREAL, XIM ,M,F0 
CALL RECORD (XREAL, IMac’ M,ISCALE) 
END. 
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ANNANANNANAINN 


ANAND 


A) <a 
(N© 


30 


60 


AND ANN 


AAD 


KKK UK KK AKA RAKKKKKA AA KKKAR AAR RAK AKAKKKKKER KK KEKE KAKA HK 
INTP2 IS A SECOND ORDER INTERPOLATION BASED ON A 


a 
* * 
* = 
* APPROXIMATED BY CENTRAL DIFFERENCES. THE INPUT VECTOR * 
* <X¥> CONTAINS UNIFORMLY SAMPLED DATA WITH UNITY % 
* TWEEN CONSECUTIVE SAMPLES. VECTOR <Y> IS INPUT * 
& H THE NEW SAMPLE TIMES AND OUTPUT WITH ~ 
* S. THERE ARE N UNIFORM SAMPLES IN <X> * 
Pe <Y 22 * 
* * 


CHOSE THE X SAMPLE TIME CLOSES TO THE WARPED 
TIME HELD IN <Y 
DO 60 LY=1,m 
CT = 100 
DG 10 IX=1,N 
IXTIME = If - 
RTABS = ABS (DT ! 
IST = PLOA (I TIME) - Y (IY) 
DIST = ABS(DIST 
IF (DTABS .LT. DIST) GO TO 25 
Dt = Y(fY) ~ IXxTIME 
CONTINUE 
IXTIME = IXTIME - 1 
jt = -1 
DO 30 J=1 
J4 = se - 1 
CaNTENUE 
TIME =_Y(I¥) 
YEE) = CN (X3, TIME, IXTIME) 
CONTINUE 
RETURN 
END | 


FCN IS THE INTERPOLATION RULE. 
FUNCTION Seer BL) 

DIMENSICN 

TX = FLOAT (1) 


COMPUTE THE COEFFICIENTS. 


A= (X(3 — eee) (1) ) 
SUB Seog eae * 


CCMPUTE THE INTERPOLATED VALUE. 


PCN = A ® (T®*®2) + B* T+C 
RETURN 
END 
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FUNCTION YLIN MAKES A LINEAR INTERPOLATION. 


AAAN 


ION YLIN (x4,YP) 
SION Xu 


° -) ¢ TO 50 
=X& (2) +¥ Pp* ‘x613)— X¥4& (2)) 
X4 (2) ~YP*® (X4 (1) ~X4 (2) ) 


1S 1) 
© 
rt 8) &S So eg 1 Pa 


FUNCTION YLAGR COMPUTES THE LAGRANGE MULTIPLIEKS 
AND MAKES THE INTEREFOLATION FOR A CHOSEN OFFSET 
FROM THE CENTRAL SAMPLE X4(2) 


AAANAAANAAN 


CTION YLAGR (X4,YF) 
ENSION X4(4 
( 


= Teen Ye Pak a A Neer -0 
YE**2Z—1) ¥ (YP-2) / 
De eae 4) )f2 


Pp—- 
scktes P**2-1) /6. 
1424 (1) #CORXG (2) +CP 1*XY (3) +CP2* X4 (4) 


MmAAANORI 
hind ORmHS 


SEKKERKMUKKKAKKKMKEKAKKHKKKAKKHKKKKKKKKKARKKEKKHKAEKKKEKKKEKKKKKKE * 
* SUBROUTINE NUPTS CALCULATES THE M EXPONENTIAL SAMPLE * 
* ECINTS FROM THE N UNIFORM SAMPLES OF THE EXISTING * 
* SPECTRUM IN PREPARATION FOR AN INTERPOLATION. THIS #*® 
* EXPONENTIALLY SAMPLED SET OF POINTS ARE STORED IN : 


* THE INPOT VECTOR X. 
PereTirerer Pero sfc at Geely £20059 064000es00s x9 00000 00Ke es 


ANAAAANNANA 


C 0 


190 
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ere ee ee ee cars ccc. 
THIS IS A STUB PROVIDING A TEST SKIN RETURN 
. TC SIMULATE A SKIN RETURN. NORMALLY THIS SUBROUTINE : 


* WILL ACTUALLY SAMPLE A SKIN RETURN. 
Hee EH TE HA He ae Ha KA ARK KE RK REAR KK ARE REKKEE & 


AANAAND 


SUBROUTINE SAMP (XR, XI,N) 
DIMENSION XR(N) ,XI(N) 


DESIGNATE SCALE FACTOR AND UNSCALED TIME SdIFT. 
= 16. / 32. 
= ae 0 A 32e 
= 1.,0/SCALE 
Sabri 


5 + 
1 
OnE Oe 


AAN 


Meta rtn iin 


CALCULATE IHE TIME OF THE SAMPLE. 
SK= (FLOAT (I)-1.0) /FLOAT (N1) 
EUILEC THE TEST SKIN RETURN. 

S Pest) CASE A DOUBLE PERIMID. 
K~T0O) ®SCALE 


GO {TO 10 
GO TO 50 


a GO TC 


AAA AAA 


deat i) ne 


WAmMerAae hela moe 


q +W 
-. EO 
0 TSCA 
* 10.0 - (TSCALE ~ TQ) * 10.0 


q 
2 
a 
3 
ic 


Mme 


(C 
(S 
~£ 
~—E 
3 
E ei 
E .G 
E .L 
T = 


G WAVE FORM IS A SQUAKE WAVE 16 SAMPLES 
AROUND THE SIXTEENTH SAMPLE. 
eee DONE ABOVE WILL EFFECT THE 


AAAAW 
2eKO © a a ae 
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( RO NESCALE 


.LT. -EDGE) GO TO 1 
<GI. +EDGE) GO TO 1 
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00 
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ANAAANIANAANQNANRANIAN 


50 


100 


KAKMRKKAKKKRKSEKAKAKKKKARKAKEKRKKHKEKRKKKKSEKKKKKKKEKKKKKKRKKKKKK 
THIS SUBRGUTINE RECORDS A SELECTED DATA SET, 
AND SCALES THE INDEX TO 32 SAMPLES (0-31). 
THE INPUTS ARE XREAL, XIMAG, AND N (THE 
NOHB EE OF SAMPLES). £SCALE fs A CCNTROL 


* 
x 
* 
x 
AELE. * 
IF ISCALE = 0 NO SCALING IS PERFORMED * 
SCALE = 1, AXIS = 31 * 

= |] (STARTS AT T = O ) * 

< 0 SHIFTING AND SCALING OCCURS “ 
[ISC LE| = AXIS LEN. A SHIFT is INCORP- * 
> * 
x 

«x 

= 


NO SHIFT (IE., THE AXIS STARTS AT QO.) 
MRR KKM AH KKM HHKKKK HRA KRKK AKKAAKAKRKEKKKEKKKKKAHKKNHEK 


SUBROUTINE RECORD (XR,XI,N,ISCALE) 
XR(N) ,X1(N) ,REC (200) 


HPHRURRHRHHEMRHHHEH HH K H 


li tin 
— + 
ue © 
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0) GO TO 40 
Isc LE) 


0) SHIFT = 0 
N, SCALE) 


S/ (FLOAT (N) -1. 
E855 KR ( ye#2) 
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AANAANAANA 


90 


100 


ANNANAAA 


90 
100 


RK KEKRKKKKKRKKEKRKKKAKKCKAKKKKKKKKKKKKKEKRKKKKKKKRKKKKKS 
* THIS CORRECTION SUBROUTINE USES ONE OF 

* THE SIMPLER CORRECTIONS FOR THE MELLIN 

* TRANSFORM. THE CORRECTION IS A PURE IMAGINARY. 


* 
Ps 
ote 
* 
CORRECTION = -FO/OMEGA . 
3 
+ 
* 


THEN A MODIFICATION [IS MADE TO THE ENTIRE 
TRANSFORM BY A MULTIPLICATION BY OMEGA. 
RRR BKMKMKRKKKRKKK TERE RAKKAKSEERRKAARARKKEKRAKKKKKK 


SUBROUTINE CORCT1(XR,X1,N,FO) 
DIMENSICN XB(N) ,XI (N) 


DO 100 I=1,N 
OMEGA ='FLOAT(I) - 1. 


XR (I) ae ‘ OMEGA 
Til) = (XL (1) - oes OHEGA) * OMEGA 
t2 {2 “= it) / ) 

XI(I) = FO 


CCNTINUE 


RETURN 
END 


RKEAKRRKRKKKRKKKKKKKKK KKAMABK SKK KRHA KKAKKHKARKKKKKURRKRKHK HE 


* THIS CORRECTION APPLIES THE MORE COMPLECATED * 
* EXPRESSION a 
CORRECTION = PO/2 + JCOT (FO/OMEGA) * 
* AND THEN MODIFIES THE ENTIRE TRANSFORM BY 1/OMEGA, * 
RKAKRKKKRKKARKRKKEAKRECKRHKAKKKAAKKKAKAKKARKARAKKKKRKKRKEKKKEKSE KX 


SUBROUTINE CORCT2 (X&,XI,N,FQ) 
DIMENSICN X(N) ,XI (4) 


PO 100 I=1,N 
CMEGA = FLOAT(I) - 1. 


XR{I) = (XEB{I) + FO/2.) *® OMEGA 

IF(I - EQ. 1)}GO TO 90 

ET ak fxr @) ~ (F0/2.) *COTAN (OMEGA) ) OMEGA 
GO 10 

XI(I) = FO;s2. 

CONTINUE 


RETURN 
END 
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APPENDIX B 


ee ee ee ee ee ee ee eee 


Se et et 


2 


MAJOR DATA STRUCTURES: 
<HFMD> 


THIS PROGRAM, FOURIER DIRECT MELLIN, TAKES AN INPUT 
WAVEFORM FRCKH fo 
FOLLOWED BY A MDMT CUTPUTING THE FEATURES TO LOGICAL 


x 

LOGICAL DEVICE 2, PERFORMS AN FFT : 
LATER PLOTTING BY MELPLT. THE MAJOR * 

U ae AND SUBRCUTINES ARE LISTED BELOW : 
D M * 
N * 
x 


Ss 


OO! 
Ht Hit baO 


EIR APPEARANCE. THE SUBROUTINES ARE 
ORE DETAIL WHERE THEY ARE ACTUALLY 
P OGRAG. 


R 
R 
I R 
* EA a HA He AR A A OE 


- THE INPUT WAVEFORM seoet 
<XFM> - THE MELLIN TRANSFOR (ae PLEX) 
<CFHI> - REAL MELLIN COEFFICI 
<SPHI> - IMAGINARY MELLIN COEFFICIENTS 
<STAND> - AN ARRAY HELD FOR LATER COMPARISON 
OR OTHER USE. 
<PRID> - AN ARRAY NORMALLY USED TO HOLD REAL 
DATA TEMPORILY. A WORK SPACE. 
<IXT> - XK AXIS TITLE FOR PLOTTING 
<IYT> - Y AXIS TITLE FOR PLOTTING 
<KEY> - NUMBER OF PLOT THIS GRAPH 


Gee 


VE - READS AN ARRAY etn LOGICAL DEVICE 
AND FILLS ZEROS TO MAKE A TOTAL OF 2 
SAMPLES. THE OOTEUT IS IN <WFM>. 

PFT - AW FFT BLOCK 

COEF ~- COMPUTES THE MELLIN TRANSFORM SAMPLE 

GHTS. THESE ARE COMPLEX 

BERS WHOSE REAL AND IMAGINARY 

TS ARE STORED IN <CPHI> AND 

HI> RESPECTIVELY. 

S A MODIPIED DIRECT MELLIN 

NSFORM TO AN INPUT WAVEFORM 

TING THE OUTPUT IN <XFM>. 

ALGORITHM IS BASED ON A FIRST 

D DIFFERENCE. 


2 
56° 


DMTM ~ APP 
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A ONODIPTED sHELLIN TRANSFORM 
DIFFERENT THAN BUT ALSO BASED 
ON THE SECOND piFFERE ENCE. 

CDMT + APPLIES A MODIFIED MELLIN TRANSFORM 
JUST AS DMIM ABOVE, EXCEPT THAT THE 
CENTRAL DIFFERENCE IS USED. 


bd 


i= 


ipl 0 
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SIME - APPLIES A MODIFIED MELLIN TRANSFORN 
USING A BACKWARD DIFFERENCE AS IN 
DMTM, EXCEPT THAT THE INTEGRATION 
IS BY SIMESON'S RULE INSTEAD OF THE 
TRAPEZOIDAL RULE. 

XAB - TAKES THE MAGNITUDE OF THE COMPLEX 
TRANSECRM <XFM> AND PUTS THE 
MANITODE IN A SPECIFIED VECTOR. 

STOX - NORMALIZE S A VE BY I 


ITH A TITLE 
E 4. 


ANOTHER. 
LE &/OR SHIFT 
OUTPUTS TO PECIFIED ARRAY 
ONC ORDER SPLINE INTERPOLATION. 
DES TWO CLOSED FORM SOLUTIONS 
RIFYING THE MELLIN ALGORITHMS. 
THE PLOTS ON THE BASIS OF 
LLING PROGRAM. 
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SO THE MAIN PROGRAM STARTS! 

DIMENSICN P&T (256) ,WFM (256) ,STAND (256) , IF (5) 
CCMMON pet ale 2) 65 CPHI (256, 135) 5 SUHT (256 pz oe PL, 
(10f bY 

DATA IF/* FR, "EQUE',* NCY ‘,° ae t/ 

FI = 3. 141592654 


oe MANY WAVEFORMS ARE TC BE TRANSFORMED? 
onuet?i) 0 


FORMAT 

KTMS IS THE NUMBER OF TIME SAMPLES (INCLUDING 
ANY ZERO FILLING IT IS A POWER OF TWO 
FCR THE CONVENIENCE OF THE FFT. 

MPTS IS THE NUMBER OF SAMPLES INPUT TO THE 
MELLIN TRANSFORM BLOCK. THE COEFFICIENTS 
ARE COMPUTED NOW. 

NU = & 

NIMS = 2#*NU 

MPTS = NTMS/2 

FCRMAT (14 

NCOZP = NTMS 

CALL COEF (NCOZF,MPTS) 

SET UP THE LOOP FOR THE NUMBER OF WAVEFORMS 
TC BE PROCESSED. 

DO 500 IWAVE=1, NUMWEM 


GET THE NEXT INPUT WAVEFORM. 
CALL WAVE (WEM,NTMS) 


ZERO THE <STAND> VECTOR TO BE USED AS THE 
IMAGINARY PART OF THE NIMS TIME SAMPLES. 


CO 100 Be Oar 
STAND ap = 
CONTINU 


CALL STCW(WEM,NTMS) 

TAKE THE FFT. 

CALL TLTLE (TF) 

CALL FFT (WFM, STAND, NTMS,NG) 


TAKE THE MAGNITUDE AND EUT IT INTO THE COMMON 
WAVEFORM <WEM> 
DO 200 I=1,NTMS 
WEM(I) = WEM Bitpen2} + STAND (I) **2 
= SCBT (WFM (I)) 


F 

ONT INU 

ALL STOW (WFM,NTMS) 
AK 

S 


4 


E THE MELLIN TRANSFORM OF THIS SPECTRUM 
ING THE FIRST HALE OF THE FFT SAMPLES, 
OTHERWISE KNOWN AS MPTS=NTMS/2. THESE ARE THE 
CNLY UNIQUE VALUES. 
CALL DMTH (#4 -MPTS, NCOEF) 
CALL XAB(PRT EY 
CALL STC (PEE NCOE ) 


jan? 
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AND 
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MWANAAND 


CALL 


USING 
SIMPSON 


L THE SECONT ORDER 
BOTH COMPUTE MELLIN 

EB SECOND DIFFERENCE APPROXIMATION 

OF THE FIRST DIFFERENCE APPROXIMATION 

- OTHERWISE THE APFROACH IS THE SAME. 


a b> 8] 

mC 

Gir 
te 


KAB (PRT, NCOEF) 
STCW (PEL ,NCOEF) 


SMT2 (WFM,MPTS, NCOEF) 
XAB(PRT, NCOEP 
STCW (PRT, NCOEF) 


CDMT WHICH USES THE CENTRAL DIIFERENCE 
FOR APFROXIMATING THE TRANSFORM. 


CDMT (WFM,MPTS, NCOEF) 
XAB TB RE NCOEE) 


giz (BE: Reoee) 


STC NCOE 


SIMP WHICH COMFUTES THE MELLIN TRANSFORM 
THE FIRST DIFFERENCE ALGORITHM AND 
*‘S RULE TO CCMPUTE THE MODIFIED 


MELLIN TRANSFORA. 


CALL 
CALL 
CALL 


END THE 
OUTPUT TO LOGICAL DEVICE 3 
TITLE INFORMATION PROVIDED BY 
FOR PLOTTING 


STAY 


CONTI 


CALL 
CALL 


STOP 
END 


SIMP (WEH, APTS, NCOEF) 
XAB(PRT BF) 
STO (pate COE ) 


TRANSFORM LOOP. THE TRANSFORM HAS BEEN 

AND PREPARED WITH 

LOGICAL DEVICE 2 
WITH PROGRAM MELPLT FORTRAN. 

IN THE LOOP Ir MORE WAVEFORMS ARE AVAILABLE. 


NUE 


CFORM (PRT, NCOEF) 
CFORM (ERT, NCOEF) 


inks 
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* THE COEF SUBROUTINE COMPUTES THE MELLIN 

* COEFFICIENTS IN TO COMMCN ARRAYS CPHI AND 
* SPHI THAT REPRESENT THE REAL AND IMAGINARY 

* PARTS RESPECTIVELY. THE TERMS ARE COMPUTED 

* BY THE FORMULA: 

* PHI(I,J) = J**S , WHERE S = A NORMALIZED 

* ESCRETE RADIAN FPREQUENCY. 

* We He Ae Ae Me he ae De ae ac he ake IC HE ae A ae Ke 9 
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SUBROUTINE COEF (NCOEF, NETS) 
COMMON KEM (256, ) CEH (256,128) , SPHI (256,128) ,PI, 
VEXT (10) LY (10) BY 
ro 100 £ = 1,NCOEF 
RI = FLOAT (Tf 
CMEGA = 2. PI * RI / 36. 
DO 20C J = 3,NPTS 
RJ = FLOAT (Jj 
CPHI (I,J) = COS(OMEGA * ALOG(RJ)) 
SEAT {T,J) = SIN(OMEGA * ALOG (RJ) ) 
200 contiNnbE 
100 CONTINUE 
RETURN 
END 
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* THE DMTM SUBROUTINE PERFORMS A DESCRETE MELLIN * 
* TRANSFORM ON THE ARRAY a ee eae * 
* FOR ONE MELLIN FREQUENC « 
x XFM (I) =SUM(K=1 TO NETS) (WEM (121) -WFM(Z) ) #Ke#S : 
* THE COMPLEX COMPONENTS CF K®*S ARE COMPUTED PRIOR * 
* TC CALLING DMTM AND STORED IN THE COMMON ARRAYS * 
* CPHI AND SPHI AS REAL AND IMAGINARY PARTS * 
* RESPECTIVELY. THE ALGORITHM IS BASED O * 
* FIRST DIFFERENCE. THE TRAPEZOIDAL RULE IS USED * 
* FOR THE INTEGRATION. THE COMPLEX OUTPUT FOR * 
* THE TRANSFORM IN IN THE COMMON ARRAY <XFM>. “ 
EEKKKKKHEKKKKEKSEKAKKAARKEKKEKKKAKKSKGAHKKKKAKKKHAKKKKKKKKKKKKES 

SUBROUTINE DMTM(SAMP,NPTS,NCOEF) 

DIMENSION FS care ) 

COMMON XPM (256,2) ,CFHI (256, 128) ,SPHI (256,128) ,PI, 

Aon Oy oon’, LLIN*,* FRES,*QOUEN*,*CY */ 

CALL TITLE (ID) 

DC 100 I=1,NCOEF 

ten tr’ =0.0 

XEM (1,2) =0.0 

N1 = fpts - 1 

DO 200 J=1,N1 

J1 = J + 1 

XFM(I,1) = CPHI (I,J) * (SAM P tS , ~- SAMP 3143 + 222 (1, 

XFM(I,2) = SPHI(I,J) * {SAM - SAMP{J1)) + XFM(I, 

CONTINUE 

CCNTINUE 

RETURN 

END 
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200 
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* SMT IS A SUBROUTINE THAT COMPUTES A NUMERICAL 
APPROXIMATICN TO THE MELLIN TRANSFORM AS DOES 
CMTM ABOVE. SMT USES A TRAPAZOIDAL APFROXTMATION 


7 
t 


C COMPUTE THE TRANSFORA, B3UT USES THE SAME 


2) 
© 


z 
See A MATRIXES <CPHI> AND <SPHI> CONTAINED 
TrTTTrTTerrr fe err errerert error irirrrere rere ere re: Shy 


SUBROUTINE SMT (SAMP,NPTS,NCOEP) 
Bean ee ToEe OMe 


10 
DATA thst rg Ete ere, Tea’, OO EN ,*CY '/ 
CALL TITTLE (ID) 


INITIALIZE THE INPUT ARRAY AND COMPUTE 
THE LCOP CONSTANTS. 

N1 = NPTS - 1 

N2= NI - 1 

SET UP THE TRANSFORM LOCE. THE OUTER LOOP 
SETS UP THE COEFFIECIENTS WHILE THE INNER 
LOOP COMPUTES THE SU WHICH ARE THE 
CCEFFICIENTS. 

DO 200 J = 1,NCOEF 

XFM(J,1) = 0.0 

XFM (3.2) = 0.0 

pe 106 I = 1,N2 

I0 =f 

Mm =I+ 1 

= tT t 2 

DELTA = SAME (10) ~ 2.* SAMP(I1) + aeifei? 
XPM(J,1) = XPM Jo 3) + DELTASCPHI Jj, 1)* 
XFM(J,2) = XFM(J,2) + DELTA*SPHI (J,1)* 
CCNTIAU 

XFM 301) = XPM(J,1) / (FLOAT (J PI/18 8.) 
XEM(J,2) = XPM(J,2) 7 (FLOAT (J) *PI/18 
XFM(J,1) = XPM(J,1) / SQRT(1+(FLOAT (J) *PI/18 
XFM (J 2) = XPM(J,2) / SCRT(1+ (FLOAT (J) *PI/18 
CONTIAG 

RETURN 

END 
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200 
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Bette HR MH 


SMT2 IS A SUBROUTINE THAT COMPUTES A NUMERICAL 
APPROKXIMATICN TO THE MELLIN TRANSFORM AS DOES 
CMTM ABOVE. SMT USES A TRAPAZOIDAL APPROXIMATION 
TC COMPUTE THE TRANSFORMS, BUT USES THE SAME 

SPO OeGn MATRIXES <CEAID> AND <SPHI> CONTAINED 
4 


KACKRARKKKKKKCKKEKKKKAKAKKRKKAAKKKAKKSARKAKKSEKKKKAKKKK KK 


S$, NCCEF) 
$3) 128) , SPHI (256,128) ,PI, 
N FR',*BQUE',*NCY '/ 


INITIALIZE THE INPUT AgaAY AND COMPUTE 


mt 
ti 
ct 
O 
© 
td 
2) 
© 
oe, 
WN 
| 
bo 
Zz, 
a | 
W— 
3 


MHRHHOe 28 CO HNCNM we 
= FT 
1 
amd 


RANSPORM LOCF. THE OUTER LOOP 
Ste cou uae WHILE THE INNER 


S 
4 


E SUM WHICH ARE THE 


re ts tH Orir 
ry 
QOr 
++ HO HE 
td Ord oo 
hou th th acm 
heim 


COaQae 
\) at 


sa + SAMP(I2)) 


- 2.% 
TA + ri ee ee, 
XFS + DFE 
ZFM + DEL 


SEI OHHAOK KO OM See 
AAI OGMAO OOnm tw 


co 4 onl a 
OH 
lt ttt Hrom 


COLOR 


ae 4 Ss 
XPM (J,2) / S 


Meo 
mimi © 
=e = 
out bd 


QRT PLOAT (J ** 2 
eaer PLOAT (3 SeBE 1g.) #*) 


FOr CC, <p HCL Oi 
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AMD 


100 


200 


RHEKAKKKCKEKKKKKKRKKAKKKKAMKKAAKAKKKKKKRKARKERKKKARKEKKKKKAKKEKKE 


HRD EHH XH 


THE COMT SUPROUTINE FERFORMS A DESCRETE MELLIN 
TRANSFORM GN THE COMMON ARRAY WEM. THE FORKULA 
FOR CNE MELLIN Fa ose VALUE I53 

XFM(X)=SUM(K=1 TO NPTS) (WFM (I +1)-WFM (I-1) ) *K¥*S 


* 
ie 
* 
* 
* 
‘HE COMPLEX COMPONENTS OF K**S ARE COMPUTED PRIOR * 
SPHI AS REAL AND IMAGINARY PARTS ” 

x 


Y. THE CENTRAL DIFFERENCE IS USED. 


eTTrTT TET TTT TTrrTTTeTTTrcrrr rrrTTtt Tr Teee TT. f 


SUBROUTINE CDMT(H, Rees NCOEF) 


DIMENSION #F 


COMMON XFM 
1IXT (19 (txt 


DAT 


CALL TITLE (1D) 


NPTS D (5) 
38 rales CPHI (256,128) ,SPHI (256,128) ,PI, 


ne), ELIN', " pee’ ,*QUENY ,*CY '/ 


SET UP DERIVATIVE OF INPUT VECTOR <H(I) > 
IDENTIFIED HERE AS <G>. 


NG = NPTS - 1 


DO 200 J 
XFM (J,1 
XFM(J,2 


OMEGA = 


ce 100 I 
ImM1 = I 
IPt =I 


Huu 


= 
=e 


+ 2 


1, NCOEF 


FLOAT (J) * PI / 18. 


1,NG 


DH = (H(IP1) - H(IM1)) / 2. 
COMPUTE THE J-TH COEFFICIENTS BY THE SUM. 


XPM (J,1) 
XFM (J, 2) 


CCUTINUE 
CCLOK = 
AFM (J,1 
ATA (J, a2 
CONTINUE 


BETURN 
END 


lb lie 


ge ,1) + DH * CPHI(J,I) 
J,2) + DH * SPHI (J,L) 


XFM(J,1) * COLOR 
XPM vip 2) * COLOR 
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fe ok 
IS SUBROUTINE USES SIMPSON'S RULE TO APPROXIMATE * 
‘HE MODIFIED MELLIN TRANSFORM. THE MODIFICATION si 
ane UENCY TIMES THE FREQUENCY DERIVATIVE OF - 
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rer. 
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SUBROUTINE Bee ibe d Ets NCOEF) 
CIMENSION H ( PTS eID (5) 
COMMON XF i 38 2 ,CPHI (256, 128) , SPHI (256, 128) , PI, 
VERE (19) EE 10) (KEY 

Et", *LLIN',* FRE',*QUEN',*CY '/ 
CALL ITTLE (ID) 


wets -1 
weIsS = 2 


DO 100 leat 


= 2.0 
CMEGA = PI * FLOAT(J) / 18. 
CCLOR = 1 


CC 200 I=1,N1 

IM1 =I 

10 = 1+f1 

IO ae ~-GT. 3.) GOTC 67 


> 
rx{ 
ie 
C, 
% 
— 
it te 


feaTey) 

wn) 
Fx 
YN 
+ 
tx 
hy 

tt 

N) 
® 


200 CONTINUE 


100 CONTINUE 
SETURN 
END 


TZ 
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* THE XAB SUEROUTING TAKES THE MAGNITUDE OF THE * 
* COMMON COMFLEX ARRAY XFM(256,2) AND PLACES * 
; seine VALUES IN OUTPUT VECTOR <XHiG> FOR LATER . 


NG 
wT eTTTTTEoT TT Tt Tt TTT ETS ST eT et CCST SCLC LT ST ST EL eT ett 


ANNNAAINNAN 


SUBROUTINE XAB XHAG ,NPTS) 

DIMENSION XMHAG (NPTS) 

CCMMON FH (258, 2) CPHI (256, 128) , SPHI (256, 128) ,PI, 
EXT (10) gt 6Y 


{ 
XMAG (I) = SQRT(XFM (I, 1) #¥2+XPM(L, zp **2 
100 CONTUAUE Seed aU, ee NE) 
RETURN 
END 


22 





G 
C SFKAREKBKKKKKKKKKKKKCKRKEKKSAKREKAKKKARKEKKSCHKAKKKKKKKARKKKKKKKX 
C * THE STOW SUBROUTINE STORES THE INPUT ARRAY PRI(NPTS) * 
C ™* INTO THE LOGICAL DEVICE 2 FOR LATER USE IN PLOTTING. * 
C #® NPLOT NUMBERS THE ELOTS FOR LATER IDENTIFICETION * 
C * AND A PLOT TITLE FOR THE HELLIN FREQUENCY IS ADDEb- * 
C ™* FOR CONVENIENCE. * 
Go + THE MODULUS OF THE * 
C * TRANSFORM TAKEN SCALED {70 UNIT MAGNITUDE, AND me 
C¢ * CUTPUT TO LOGIC CAL * 
a * INEUT: - UTS. BE SCALED TO 1 AND WRITTEN * 
a TO LOGICAL UNIT 2. * 
a = NPLOT - THE NUMBER OF THE PLOT x 
eG KEY NUMBER OF CURVE THIS PLOT * 
Cc Rae AE AOS bh eck KA ERAAE RAR AREER SES ARES ROR ERR RKAES 
c 


SUBROUTINE STOW (PET, NPTS) 

DIMENSICN ERT {NET i 

COMMON XE 4 (25 2) CPHI (256, 128) , SPHI (256,128) ,PI, 
1IXT (10) ,LYT (iC) ,K 


WRITE (3,13) NPTS 
is 3,13 KEY 


CKK E AREY + A)GQ TO 12 
C GRIVE Axis LAB 
WRITE (3,10) (XT (T) ,1=1, 10 
WRITE (3-10)-(1Y@ eta] *10 
10 FORMAT (104). 
12 NELOT = NELOT + 1 
BPET = 0.0 
13 FORMAT (IG) 
DO 100 I=1,NPTS 
Iz {Pat (I) °.GT. BPR1) BERT = PRT (I) 
100 CCNTING 
IF (BPRT -LT. 00001) BPRT = 1. 
DO 200 I=1,NPTS 
FET (Z) = PRT(I) / BPRT 
T= ELOAT (T) 
WRITE (3,20) T, BRI (I) 
20 EQRMAT (2 (3 ,F9.5)) 
200 CCNTINUE 
RETURN 
END 
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* THE HOLD SUBROUTINE TAKES THE INPUT FILE - 
* <FILIND> AND STORES IT IN THE OUTPUT FILE * 
* <FILOUT> FOR TEMPORARY STORAGE. THE FILE * 
* <FILIN> REMAINS UNCHANGEL, * 


4 HH NA He ee Re i a eR a a A A a CK Ce 
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SUBROUTINE HOLD (FILIN,FILOUT,NPTS 
CIMENSION FPILIN (NPTS),FILOUT (NPTS 


OC I = 1,NPTS 
mL oer At) = FILIN (I) 
100 eee E 


EREKKARKKAKKKKEKKKEKKKKKKKAKKAKKKAKE KCKAKKKKKKKKEKKKAEKKEKKE 
* THE ALTER SUBROUTINE WILL ALTER THE COMMON ARRAY * 
* <WFM> AS SPECIFIED BY THE INPUT VARIABLES * 
* <SCALE> AND <SHIFT>. THE ALTERED WEFM IS OUTPUT : 


* IN THE VECTOR <ALT>. 
OTT TILT LL TTL tT TTP ttt ett tt etree ett et ee et 


AAAQAQNAANN 


SUBROUTINE ALTER (ALI, ¥FM,SCALE, SHIFT, NPTS) 
DIMENSION ALT (NBTS) gHEM ( {NPTS) TOLD 256) TNEW (256) 
8» 2) CE HI (256, 128) , SPHI (256,128) ,Pl, 


erat 
MOAOORAOHOA 
2 tt O2OObKO 
OAMamr He 


woo? 


=1 pueTs T (T) 
= TOLD(1) / SCALE + SHIFT 


E 
NIP3(WFM,TOLD,ALT,TNEW,NPTS) 


RH Oe O 
SHoHBZHHOoOose 
oeo~ 


100 


GSt*4 m0 & 


3 @] 
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ecm 


INTP3 IS A SECOND ORDER INTERPOLATION BASED ON A 
ECLINOMIAL FOUATING TO FIRST AND SECOND DERIVATIVES 
APPROXIMATEL BY CENTRAL DIFFERENCES. THE LNPUT 
VECTOR <XO> HAS OLD SAMELES AT TIMES IN <TO>. THE 
NEW SAMPLE TIMES ARE INPUT THROUGH ARRAY <TN> AND 
THE COMEUTED SAMPLES AT THESE TIMES ARE OUTPUT LN 
ARE OUTPUT IN THE VECTOR <XN> AND <NFM>. 
TT ETETTTOTELET TLE ELITE TL TTT rrTrrTrrre, 
SUBROUTINE TTP 3 (x0 TO, XN,TN,NPTS) 
DIMENSION £313) O Apts) fo (Apts) XN (NETS) , TN(NPTS) 
COMMON XFM(256,2) ,CFHI (256,128) , SPHI(256,128),PI, 
IXT (10) ,LYT (10) , KEX 
CHOSE THE <TO> SAMPLE TIME CLOSES TO THE WARPED 
TIME HELD IN <TN> 
pC 60 I=1,NPTS 
DT = 100 
XPTS = FLOAT (NPTS 
LE ( (TN (I) »AND. (TN(I) .LE. XPTS)) GO TO 5 
xn(t) = 0.0 
GO TO 60 
ro 10 J=1,NPTS 
DTABS = ABS{DT) 
LIST = TO(J). - TN (ZI) 
DIST = ABS(DIS *} 
IF (DTABS “LT IST) GO TO 25 
DI = IN¢I) - 
CCNTINUGE 
JTIME = 3 - 1 
Ji = -1 
DO 30 J=1,3 
J1 = J + JTIME-1 
IF 4 (31 .GE. 1) .AND. (J1 .LE. NPTS))GO TO 29 
X3(3) = 0.0 
gS ay a (31) 
ceNtinoe 
TIMN = TN (I). 
TIMO = TO JIIME) 
XN (I) = PCN (x3, TIMO,TIMN) 
CONTINUE 
[LO 500 I=1,NPTS 
CONTINUE 
BETURN 
END 
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ANN ANAND 


ANNA 


FCN IS THE INTERPOLATION RULE. 


FUNCTION FCN(X,TO,TN) 
DIMENSION X (3) 


COMPUTE THE COEFFICZENTS. 


A = = 28 x (2), - o*t ) 
B= x 3} > EEL hor? 208 A *% 70 


COMPUTE THE INTERPOLATED VALUE. 
FCN = A *® (IN**®2) + 3 ® TN + C 


KETURN 
END 
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@ KKRAKBRKRKKKKKAKKRKKCKKAKAKRKERAAEKAKARAK KCKERAKKAKKHK EKRKKKKKKAKKKEE 
a * THIS SUBROUTINE IS THE RESULT OF A CLOSED * 
cS * FORM CALCULATION OF THE CONTINUOUS SINC#¥*2 * 
iS * EFING TAKEN FROM THE TIME DOMAIN THROUGH THE * 
& * ENTIRE FOUBIER-MELLIN PROCESSING. THE ALGORITHA * 
c * IS USED TO VERIFY THE FM PROCESS. THE OUTPUT * 
c * FEATURE SPACE SHOULD BE IDENTICLE TO THAT PRODUCED 7 
c * BY ANY OF THE MELLIN ALGORITHMS. FOR THIS REASON 7 
2 * THE SAMPLE POINTS ARE SYNCHRONIZED WITH THOSE USED os 
C * ABOVE. <NCOEF> FM COEFFICIENTS ARE USED. THE * 
C * CCMMON VARIABLE <KEY> MUST BE SET TO 99 PRIOR TO = 
C * CALLING TQ GET T SINC*®*2 OUTPUT FEATURE SPACZz. * 
C * TQ OUTPUT THE CLOSED FORM SOLUTION FOR A RAMP IN * 
c * THE PREQUENCY DOMAIN, <KEY> MUST EQUAL 100. * 
e shee 3H a he 2h eh ae fs 3 he fe fe ae he a hehe He He He He He hee 3 hee He he a oe 2 he 2 ae ee i ae A a ae ae a a a a 
S 

c 


SUBROUTINE CFORM(CF,NCOEF) 
DIMENSICN CF (NCOEF) 


CCMMON XFM (256, 2) , CPHI (256,128) , SPHI (256,128) ,PI, 
11XT (10) ,LYT (10) , Kz 


IF(KEY .EQ. 100) GO TO 200 
LP(KEY .NE.99) RETURN 
CO 100 M= COEF 


. 
= 


XM = FLOAT (M) - 1.0 

CEiM) | "SOR raMAG)~ “anag* 
100 CONTINGE / 

GC TO 101 
200 DO 202 I = 1,NCOEF 

XI = FLOAT(I) -1. 

OMEGA = PI * x 18. 

2./SQRT (4. + OMEGA**2) 


202 CO ONT INTE 
101 CONTINUE 


READ (2,20) KEY 

READ (2, 10 IXxT T=1 
READ(2, {itt r= 
PoRnAT( 104d 
FORMAT (14) 


CALL STOW (CF,NCOEP) 


RETURN 
END 
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mM TLL sPiTmes THE Y AXIS FOR THE: PLOTTING ¥ 
* FROGEAM ACCORDING TO THE ALGORITHM JSED TO * 
* GENERATE THE FEATURE SPA * 


ANNAAY 


STITT TT TiTtTrtt tT T Tt tT Tt TTT TTrrrrrererlTTT Te TTT eTy te 


SUBROUTINE TITLE (ID) 


5 
ON XFM ( 8 2) «¢ CPHI (256, 128) , SPHI (256,128) ,PI, 
9) TYT (19) 


1 
00 
thor 2 16 (2) 
RN 


=" 


100 


128 
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